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PREFACE 


The  practical  implications  of  linear  sequential  filtez'ing  (or 
Kalman-Ducy  filtering)  theory  were  quickly  recognized  by  the  engineering 
community  as  an  important  contribution  to  -■'eal-time  data  processing.  Some 
of  its  numerous  successful  applications  have  been  made  to  aerospace  en¬ 
gineering  system.  The  typical  examples  are  orbit  determination  and  trajec¬ 
tory  estimation  problems.  Since  these  problems  generally  are  concerned  with 
continuous  nonlinear  dynamic  systems  and  discrete  observations,  nonlinear 
filtering  has  been  a  theme  of  interest  in  the  field  of  orbit  determination. 

It  has  been  demonstrated  that  the  optimal  nonlinear  filter. requires 
the  computation  of  an  infinite  number  of  moments  and  generally  its  implementa¬ 
tion  is  not  practical.  Tnis  leads  one  to  seek  an  approximate  solution  to  the 
optimal  nonlinear  filtering  problem.  Several  approximate  nonlinear  filters 
have  been  proposed  previously  and,  for  the  most  parti  these  can  be  classified 
as  one  of  two  basic  types  of  second  order  filters.  The  first  is  the  truncated 
second  order  filter  which  utilizes  a  Taylor  series  expansion  of  the  dynamic 
system  and  the  state-observation  relationships ,  followed  by  a  truncation  of 
the  third  and  higher  order  moments.  The  other  is  the  Gaussif.n  second  order 
filter  which  employs  a  Taylor  series  expansion  and  approximations  of  the 
fourth  order  moments  in  terms  of  the  second  order  moments ,  under  the  assump¬ 
tion  that  the  conditional  density  function  is  Gaussian.  The  unique  feature 
of  botli  filters  is  found  in  the  fact  that  a  random  forcing  term  occurs  in  the 
covariance  equation.  The  random  forcing  term  which  enters  into  the  covariance 
equation  in  a  linear  manner  is  considered  to  have  potential  for  causing  the 
conditional  covariance  matrix  to  be  negative  definite  over  some  non-zero 
time  interval.  This  term  is  often  neglected  in  the  modified  Gaussian  or 


Hi 


tr-unrated  second  order  filters. 


This  study  is  concerned  with  the  development  of  an  approximate  non¬ 
linear  filter  using  the  Martingale  theory  and  approoriate  smoothing  proper¬ 
ties.  Both  the  first  order  ar-u  the  second  order  moments  are  estimated.  The 
filter,  which  is  developed,  can  be  classified  as  a  modified  Gaussian  second 
order  filter.  Its  performance  is  evaluated  in  a  simulated  study  of  the  prob¬ 
lem  of  estimating  trie  state  of  an  interplanetary  space  ''ehicle  during  both  a 
simulated  Jupiter  fly-by  and  a  simulated  Jupiter  orbiter  mission.  In  addi¬ 
tion  to  the  modi!  ied  Gaussian  second  order-  filter,  the  modified  truncated 
second  order  filter  is  evaluated  also  in  the  simulated  study.  Results  ob¬ 
tained  with  each  of  these  filters  are  compared  with  numerical  results  obtained 
with  the  extended  Kalman  filter  and  the  performance  of  each  filter  is  deter¬ 
mined  by  comparison  with  the  actual  estimation  errors.  The  simulations  are 
designed  to  determine  the  effects  of  the  second  order  terms  in  the  dynamic 
state  relations,  the  observation-state  relations  and  5n  the  Kalman  gain  com¬ 
pensation  term.  The  result  of  an  extensive  simulation  shows  that  the  Kalman 
gain  compensated  filter  which  includes  only  the  Kalman  gain  compensation  term 
Is  superior  to  all  of  the  other  filters. 
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CHAPTER  1 


INTRODUCTION 

1 . 1  Preliminary  Remarks 

In  the  field  of  space  tracking  and  guidance,  one  basic  requirement 
for  spacecraft  guidance  is  the  capability  to  obtain  and  to  rapidly  process 
observations  to  determine  an  estimate  of  the  spacecraft  trajectory.  This  re¬ 
quirement  initiates  the  search  for  mathematical  techniques  which  are  compu¬ 
tationally  efficient,  but  which  possess  a  high  degree  of  accuracy.  Follow¬ 
ing  the  precise  formulations  of  the  linear  sequential  estimation  theory  by 
Kalman  and  Bucy  (6,7),  the  practical  implications  of  the  theory  were  recog¬ 
nized  and  numerous  successful  applications  have  been  made  in  the  field  of 
orbit  determination  and  trajectory  estimation  problems,  however,  these  ap¬ 
plications  generally  are  concerned  with  nonlinear  continuous  dynamic  systems 
and  nonlinear  state-observation  relationships  and,  hence,  the  linear  esti¬ 
mation  theory  cannot  be  applied  directly.  As  a  matter  of  fact,  it  is  not  a 
simple  task  to  apply  the  linear  estimation  theory  to  orbit  determination 
problems.  Usually  nonlinear  dynamic  systems  and  state-observation  relation¬ 
ships  are  linearized  about  a  nominal  (or  reference)  trajectory  under  the  as¬ 
sumption  that  the  true  trajectory  is  sufficiently  close  to  the  reference  tra¬ 
jectory,  and  then  the  linear  estimation  theory  is  applied  to  the  linearized 
systems.  Conceptually,  there  are  two  ways  to  carry  out  the  linearizations 
and  the  resulting  filters  are  somewhat  different  from  each  other.  The  dis¬ 
tinction  is  how  the  nominal  trajectory  is  chosen.  If  a  prescribed  trajectory 
is  chosen  as  a  nominal  trajectory,  *.he  original  Kalman-Bucy  linear  filter 
can  be  directly  applied  to  the  linearized  system  which  governs  the  state  and 
the  observation  deviations  from  the  prescribed  nominal  values  of  the  state 
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and  observation.  Although  this  approach  appears  to  be  conceptually  simple, 
it  suffers  from  several  points.  First,  if  the  nominal  trajectory  is  not 
close  enough  to  the  true  trajectory,  the  basic  assumption  used  in  the  lin¬ 
earization  procedure  is  violated  and  the  estimate  of  'the  deviation  from  the 
nominal  trajectory  filter  can  lead  to  inaccurate  results  and  often  diverges. 
Furthermore,  it  is  intuitively  more  appealing  to  take  the  current  estimate, 
rather  than  a  prescribed  trajectory,  as  a  nominal  trajectory  and  conduct 
the  linearization  about  the  current  estimate  instead  of  a  prescribed  trajec¬ 
tory.  In  this  case,  the  linearized  system  will  involve  deviation  in  the 
state  and  the  observation  from  the  current  estinates  of  the  state  and. obser¬ 
vation  instead  of  values  related  to  u  prescribed  nominal.  The  original 
Kalman-Bucy  linear  filter  can  be  applied  to  the  above  linearized  system.  The 
advantages  of  using  the  current  estimate  as  a  nominal  trajectory  are  that  a 
nominal  which  is  closer  to  the  true  trajectory  can  be  used  and  that  the  fil¬ 
tering  procedures  can  be  simplified  due  to  the  fact  that  all  the  propagated 
state  deviations  will  become  identically  zero.  This  concept  will  be  clearly 
discussed  in  Section  1.3.  In  order  to  represent  this  situation,  "the  extend¬ 
ed  Kalman  (EK)  filter"  proposed  by  Jazwinski  (2)  and  distinguished  from  the 
prescribed  nominal  trajectory  filter  will  be  used. 

It  is  well  known  that  all  of  the  information  about  the  state  pro¬ 
vided  by  the  measurements  is  contained  in  the  probability  density  function 
of  the  state  conditioned  upon  the  entire  past  history  of  measurements.  From 
this  conditional  probability  density  function,  one  can,  in  principle,  deter¬ 
mine  the  optimal  filter.  In  general,  the  optimal  filter  is  expressed  in 
terms  of  the  moments  of  the  conditional  probability  density  function.  Hence, 
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this  conditional  density  function  becomes  a  prime  ingredient  for  studies  of 
optimal  filtering. 

Several  authors  have  considered  the  problem  of  deriving  a  dynami¬ 
cal  equation  for  the  conditional  density  function  when  the  dynamic  state 
noise  and  observation  noise  are  both  jointly  Gaussian  and  white.  The  most 
recent  pattern  of  research  in  this  field  appears  to  have  been  initiated  by 
Stratononich  (12)  and  Kushner  (0,9).  The  formal  chnracte.  of  this  initial 
work  stimulated  numerous  studies  of  nonlinear  filtering  which  have  attempted 
to  extend,  and  to  obtain  a  more  rigorous  verification  of  these  initial  re¬ 
sults.  The  method  used  by  Strator.on  ich  and  Kushner  is  based  on  a  discrete 
time  model,  and  an  iterative  application  of  Bayes'  rule  is  used  to  obtain  a 
representation  of  the  conditional  density  function.  The  solution  of  the 
continuous  time  problem  is  obtained  by  a  limiting  process. 

Although  the  central  ideas  and  methods  were  all  supplied  by 
Stratonovich  and  Kushner,  and  most  other  papers  in  this  area  are  just  con¬ 
cerned  with  extentior.  of  these  basic  ideas,  Bucy's  ('1,5)  approach  to  the  op¬ 
timal  nonlinear  filtering  problem  is  rather  unique  and  more  mathematical  than 
Stratonovich  and  Kushner's.  However,  the  results  which  he  obtained  for  the 
Gaussian  state  and  observation  noise  case  were  identical  to  those  of  Kushner. 
An  important  intermediate  result  of  Bucy's  work  is  that  of  a  representation 
theorem  which  demonstrates  how  the  posterior  cunditionaJ  density  function  at 
some  instant  of  time  can  be  represented  as  a  function  of  the  apriori  density, 
P(X(to))  and  the  conditional  expectation  of  an  exponential  function  of  the 
observational  data  over  the  time  interval  (t  ,t)  . 

In  addition  to  the  above  research,  Mortenson  (1"),  Cox  (25)  and 
IVtrhmondy  et  al.  (21)  approach  the  nonlinear  filtering  problem  from  the 


control  theoretical  point  of  view.  They  formed  a  likely  hood  function  and 
maximized  the  function  in  various  ways,  e.g.,  dynamic  programming  or 
Pontryagin's  maximum  principle. 

Fisher  (23)  presented  a  unified  and  compact  development  of  the  non¬ 
linear  filtering  problem  for  a  broad  class  of  Markov  signal  processes,  by 
making  use  of  the  characteristic  function  technique.  The  idea  of  approaching 
the  continuous  time  nonlinear  estimation  problem  from  the  Innovation  process 
approach  was  suggested  by  Frost  (33). 

There  have  been  a  number  of  associated  approximation  methods  de¬ 
veloped.  Noteable  among  those  are  those  of  Bucy  (4,5),  Kushner  (10), 
Jazwinski  (1,2,3),  Bass  et  al.  (18,20),  Schwartz  (19)  and  Athans  et  al.  (27). 
Most  of  the  references  cited  above  utilize  techniques  that  are  closely  re¬ 
lated  to  the  methods  introduced  by  Kushner  (10)  and  Bucy  (4),  namely,  Taylor 
series  expansion  and  the  assumption  of  a  Gaussian  density  function  or  Taylor 
series  expansion  and  truncation.  Utilizing  the  Taylor  series  expansion  tech¬ 
nique,  there  are  two  basic  types  of  second  order  filters  which  have  been  de¬ 
veloped.  Tirst,  it  is  assumed  that  the  third  and  higher  order  moments  are 
negligible.  The  resulting  filter,  referred  to  as  the  truncated  second  order 
filter,  was  developed  by  Jazwinski  (1,2,3)  and  independently  by  Bass,  Norum 
and  Schwartz  (20)  who  extended  the  idea  of  Bucy  (4)  to  the  arbitrary  n-dimen- 
sional  case.  Schwartz  (19),  Jazwinski  (3)  and  Fisher  (24)  independently  de¬ 
veloped  the  Gaussian  second  order  filter.  In  this  approximation,  the  fourth 
order  moments  are  approximated  in  terms  of  the  second  order  moments  under  the 
assumption  that  the  conditional  density  function  is  Gaussian. 

A  significant  feature  of  both  the  truncated  and  the  Gaussian  second 
order  filters  is  ths  presence  of  a  random  forcing  term  in  the  covariance 
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equation.  The  presence  of  the  random  forcing  term  5s,  in  principle,  justi¬ 
fiable.  However,  there  has  been  considerable  controvi  ,’sy  associated  with 
the  presence  of  this  term. 

The  term  enters  with  a  plus  sign  in  one  filter  and  with  a  minus 
sign  in  the  other.  Furthermore,  the  term  enters  in  linear  manner  and  there 
is  a  possibility  that  a  negative  variance  may  result  due  to  the  sign  of  the 
observation  residual.  These  considerations  suggest  r  ,'nsidering  a  compromise 
between  the  truncated  and  the  Gaussian  second  order  filters.  Jazwinski  (3) 
dropped  the  forcing  term  in  the  covariance  equation  for  the  compromise  and 
defined  the  modified  truncated  second  order  (MTSO)  filter  and  the  modified 
Gaussian  second  order  (MGSO)  filter,  respectively. 

Athans,  et  al.  (27)  developed  the  modified  Gaussian  second  order 
filter  using  an  assumption  based  on  an  intuition  argument  and  applied  the 
filter  to  a  simple  one-dimensional  free-fall  reentry  problem  with  range  type 
of  measurement.  The  result  of  the  simulations  indicates  considerable  promise 
for  the  MGSO  filter. 

In  this  report,  the  modified  Gaussian  second  oi’der  filter  was 
developed  rigorously  on  the  basis  of  the  Martingale  theory  and  smoothing  pro¬ 
perties  of  Loovo  (34).  The  resulting  algorithm  was  applied  to  a  study  of 
the  Jupiter  fly-by  and  the  Jupiter  orbiter  missions  using  range,  range-rate, 
star-planet  'nd  sun-planet  angle  measurements. 

1.2  Kalman-Bucy  Filter 

Consider  the  linear  dynamics  system  described  by  the  linear  vector 
stochastic  differential  equation 

dx  =  F(t)x  dt  +  G(t)da  ,  t  %  l  (1.1) 

t  t  t  -  o 


0 


Uq.  (1.1)  can  be  expressed  form, illy  as  (see  Ref.  lb) 

dx^. 


dt 


F(t)xt  +  G(t)ut  ,  t  :  t 


(1.2) 


where  x  1:5  the  nxl  state  vector,  l'(  * )  and  G(-)  are,  respectively, 
n*n  arid  non-random,  continuous  matrix  functions  of  time,  and  >  t  J 

is  an  r-vector  Brownian  motion  process  with  the  statistics 


£[dfcds[j  =  Q(t)dt 


The  r-vector  u  is  a  white  Gaussian  vector  process  which  can  be  regarded 


as  the  time  derivatives  of  . 


It  is  assumed  that  linear  observations  are  taken  at  discrete  time 
instants ,  k  : 


H(k)x]<  +  vk  ;  k  =  1,2,.. 


(1.3) 


where  y^  is  an  m-vector  of  observations,  H(*)  is  an  m'n  non-random, 
bounded  matrix  function,  and  (v^.k  =  1,2,...}  is  an  m-vector,  indepen¬ 
dent  Gaussian  sequence,  i.e.  ,  v^  --  N(0,Rk),Rk  >  0  for  all  k  .  The  dis¬ 
tribution  of  xq  is  Gaussian,  i.e.,  xq  'v  N(£q,P^)  ,  and  x  9 }  and 
{v^}  are  assumed  to  be  independent. 

The  fact  that  the  minimum  variance  estimate  is  given  by  the  condi¬ 
tional  expectation  (see  Appendix  T),  leads  to  the  requirement  that  the  con¬ 
ditional  expectation  &t/t  ~  E[x  /YtJ  for  the  above  system  bo  found.  The 
solution  to  this  problem  yields  differential  equations  of  evolution  for  the 


conditional  expectation  £  .  and  the  covariance  matrix  P  .  .  Between 

'  t+s/t  t+s/l 


observations,  these  relations  satisfy  the  differential  equations  given  in 
Cqs.  (1.4)  and  (1.5),  respectively. 
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‘I  tn/t 
(Ik 


'(c+n)*t+s/t  •  *(t)  =  Vt 


(1.4) 


d" 


=  F(t+s)?t+s/t  +  Ft+s/tFT(t+s)  +  G(t+S )Q(t+s)GT(t+s)  ,  P(t)  =  Vt/t  (1.5) 


where  the  superscript  T  indicates  th  transpose  of  matrix,  t  represents 
the  time  at  which  the  last  observation  was  made,  and  s  represents  any  time 
segment  after  t  and  before  a  new  observation  is  made.  At  the  instant  t+s  , 
namely,  immediately  after  a  new  observation  is  incorporated  at  t+s  the 
following  difference  equations  are  satisfied. 


9  - 

't+s/t+s 

^t+s/t  +  ^t+s^t+s  -t+s/t  ^ 

(1.6) 

V  = 

t+s/t+s 

Ft+s/t  "  Kt+sH(t  +  s)Ft+s/t 

(1.7) 

^t+s/t 

n(t  +  s)*t+s/t 

(1.8) 

K.  .y 

t+s 

W"T(t  +  s»H(t  +  +  »■)  +  W' 

(1.9) 

^t+s/t  *s 

the  predicted  observation  on  the  basis  of  Yt  and 

K,  , 

t+s 

is  the  Kalman  gain. 

The  solution  of  Eqs.  (1.4)  and  (1.5)  are  referred  to  as  the  apriori 
estimate  and  the  apriori  covariance  matrix,  respectively.  Meanwhile,  the 
solutions  of  the  difference  Eqs.  (1.6)  and  (1.7)  are  said  to  be  the  posterior 
estimate  and  the  posterior  covariance  matrix,  respectively.  Once  the  pos¬ 
terior  estimate  ^t+s/t(.s  and  the  posterior  covariance  matrix  Vt+s/.t+s 
are  obtained,  they  can  be  used  as  initial  conditions  for  the  differential 
Eqs.  (1.4)  and  (1.5),  respectively.  By  integrating  these  relations  forward 
until  a  now  observation  is  obtained,  the  apriori  estimate  which  is  the  op¬ 
timal  estimate  between  obsei vati.os  is  obtained.  In  order  to  initiate  this 


From  the  statistics 


procedure ,  it:  is  .necessary  to  specify  &Q/Q  and  V 

of  rondo,.,  variable  xq  ,  &q/q  and  Vq/q  are  given  as  Sq  and  Pq  ,  res¬ 
pectively. 

1.3  Linearization  and  the  Extended  Kalman  Filter. 

As  pointed  out  previously ,  the  linearization  of  the  dynamic  systeih 
aid  the  state-observation  relationships:  cannot  be  avoided,  if  the  original 
Kalman-Bucy  filter  is  applied  io  the  orbit  determination  problem.  The  de¬ 
tails  presented  in  the  subsequent  discussion  are  used  to  obtain  the  extended 
Kalman  filter. 

Suppose  that  the  equation  of  motion  is  described  by  the  following 
nonlinear  stochastic  differential  equation 

dXr  =  f(XT,T-)dT  +  G(T)dBT  ,  T  >  ro  (1.10) 

Eq.  (1,10)  can  be  expressed  formally  as 

dXT 

&T  =  f(xT,.t)  +  G(t)ut  (1.H) 

The  discrete  nonlinear  observations,  which  are  taken  at  time  instants  k  , 
can  be  expressed  as 

Yk  =  hiXk)  +  vk  »  k  =  1*2,3,...  (1.12) 

In  the  above  systems,  ST,u^,v^  and  xq  are  assumed  to  have  the  properties 
described  in  Section  1.2.  The  o-field  generated  by  the  observations  Y,  is 
denoted  by  Z t  ,  that  is  =  f ;0  <  k  <  t)  .  Substituting  XT  -  X*  t  xt  , 
and  expanding  f(XT,x)  in  Eq.  (1.10)  about  the  nominal  X*  at  each  point 
in  time  leads  to 


f (X*  +  x  )T )dt  +  G(-r)d|3 

T  T  T 


f(X*,t)dT  +  fx(X*,f  )xt<3t  +  G(;)d3T 


where  terms  involving  powers  of  x^  higher  than  the  first  >ne  are  neglected. 
For  the  nominal,  the  following  equation  must  be  satisfied 


(1.14) 

Hence,  the  state  deviation  nan  be  described  by  the  linear  time  varying 
•stochastic  differencial  equation. 


(1.13) 


dxt  =  fx(Xt’t)xtdt  +  G(t)det 


(1.15) 


The  same  procedure  can  be  applied  to  the  state-observation  relationships 
(1.12),  and  the  final  result  would  be  expressed  as  follows 


(1.16) 


\KP\  +  ’'K 


(1.17) 


for  the  observation  deviation.  Here  h^  is  the  first  partial  derivative 
of  h( • )  m.r.t.X  . 

Combining  Eqs.  (1.15)  and  (1.17)  leads  to  che  same  model  which  was 
discussed  in  the  previous  Section.  Therefore,  the  Kaiman-Bucy  linear  filter¬ 
ing  theory  can  be  applied  to  the  system  of  Eqs.  (1.15)  and  (1.17).  Hence, 
the  optimal  estimate  of  state  deviation  between  observations  is  given  by  the 


,  ^aiuSis  -hi  t&xi-i.  <-  - 
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solution  of  the  following  linear  differential  equation 

d^t+s/t 

-dr4-  ■  y*?+s>t+°>w 

At  the  observation  time  t  +  s 

*t+s/t+s  “  *t+s/t  +  Kt+s(yt+s  “  i>t+s/t) 

where 


(1.18) 


a. 19) 


^t+s/t  "  hx*Xt+sHt+s/t 


(1.20) 


Let  t  =  t  +  s  ;  then  dt  =  ds  and  substituting  these  in  Eq.  (1.14)  will 


yield 

dx* 

"dP  =  *<*?**  ♦  s) 

Combining  Eqs.  (1.18)  and  (1.21),  the  following  result  is  obtained 
d(Xt+s  +  W 

- 35 - -  *  *  s>  +  **«?„.*  ♦ 


(1.21) 


(1.22) 


Denoting  ^+s/^.  -  +  >  which  is  an  apriori  estimate  of"  X 

-  tl’l"  s 

based  on  ,  Eq.  (1.22)  can  be  approximated  as 


tts/t  _ 

ds  "  Xt+s/t  + 


(1.23) 


From  Eqs.  (1.21)  and  (1.23),  it  follows  that  the  same  differential  equation 
governs  the  nominal  trajectory  X*+g  as  well  as  the  apriori  estimate  xt 
Therefore,  selecting  the  same  initial  conditions  for  Eqs.  (1.21)  and  (1.23) 
will  lead  to  the  conclusion  that  the  no~!nal  trajectory  and  the  apriori  es¬ 
timate  are  identical.  This  situation  is  satisfied  if  the  current  optimal 
estimate  is  chosen  as  a  nominal  trajectory.  In  other  words,  if  the  nominal 
trajectory  is  updated  with  a  current  optimal  estimate,  the  apriori  estimate 


^t+s/t  °'  state  x-t+s  >  governed  by  Eq.  (1.18)  becomes  identically  zero, 
and  hence,  Eqs.  (l.ifi)  and  (1.20)  are  not  necessary.  For>  this  situation, 
the  initial  condition  for  Eq.  (1.18)  would  be  zero  and,  consequently,  the 
solution  becomes  identically  zero. 

From  the  induction  above,  it  follows  that  it  would  be  simpler  to 
linearize  the  system  about  the  current  optimal  estimate  instead  of  a  certain 
prescribed  nominal  trajectory.  This  situation  usually  occurs  in  nature. 

For  example,  when  one  deals  with.- the  motion  of  a  rigid  body,  it  is  always 
hotter  to  stick  with  the  mass  center,  which  is  the  analogy  of  the  optimal 
estimate  of  position ,  that  is ,  the  mass  center  is  nothing  but  a  conditional 
expectation  of  equivalent-point  mass.  Furthermore,  it  has  been  demonstrated, 
in  the  numerical  simulations  (30,43),  that  taking  tne  current  optimal  estimate 
as  the  nominal  trajectory  leads  to  better  convergence  characteristics  than 
using  a  certain  prescribed  nominal. 

At  observation  time  t  +  s  ,  the  optimal  estimate  ^+s/t+s  °? 

^t+s  wou-*-<*  ke  expressed  as  follows 


Y  —  v*f*  i  a 

At+s/t+s  t+s  xt+s/t+s 


(1.24) 


If  the  optimal  estimate  is  chosen  as  a  nominal  trajectory,  then 

X*  =  X.,  .  Therefore,  Eq.  (1.24)  becomes 

t+s  t+s/t  ’  1 

^t+s/t+s  ^t+s/t  +  ^t+s/t  +  \+s^t+s  ^t+s/t^ 


(1.25) 


Since  and  ^t+s/t  are  zero  ^or  case  w^ere  optimal  es¬ 


timate  is  chosen  as  a  nominal,  Eq.  (1.25)  again  becomes 


S:+s/t+s  ~  *t+s/t  +  Kt+s^t+s  +  Yt+s  "  Yt+s^ 


(1.26) 


^ri 7 ..Bj  [ 


after  adding  and  subtracting  Y*+g  -  It  can  be  easily  seen  from  Eq.  (1.16) 


Mi.it 


Y?tS  =  "<X?+s>  4  h(!W> 


t+s/t 


(1.27) 


and 


t+s 


y . ,  +  Y* 

Jt+s  t+£ 


(1.28) 


t+s/t  _ 
ds 


and 

dP 


t+s/t  _ 


Therefore,  Eq.  (1.26)  finally  becomes 
X 


t+s/t+s 


Xt+s/t  +  Kt+s^Yt+s  "  Yt+s/t^ 


(1,29) 


The  above  expression  represents  the  extended  Kalman  filter  (3)  and  can  be 

«■> 

summarized  as  follows;  between  observations,  the  apriori  estimate 
and  the  apriori  conditional  covariance  matrix  which  is  distinguished 

from  the  apriori  covariance  matrix  for  the  linear  system,  must  satisfy  the 
following  ordinary  differential  equations 

A 

dX. 


"Wt-* +  s) 


(1.30) 


=  fx(*t+s/t’t+s)Wt  '  Vs/tfi(*t+s/t’t+s)  +  G(t+s)Q(t+3)GT(t+s)  (1.31) 


respectively.  At  the  observation  rime  t  +  s  ,  the  posterior  estimate 
^t+s/t+s  anc*  t*le  Poster;i-or  conditional  covariance  matrix  ^^s/j.+s  are 
determined  by  the  following  set  of  difference  equations. 


t+s/t+s 


t+s/t+s 


+  K  (Y  —  Y  \ 

't+s/t  t+s''  t+s  t+s/t' 

(1.32) 

’  -  K  h  (X  'IP 

t+s/t  *  t+SX'At+s/-C  t+s/t 

(1.33) 
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W  ■-  h'W> 


(1.34) 


t+s 


W‘xt*twt>th*(W>Whx(*t«/t> -  ■W'1  «-35> 


or 


K 


t+s 


^t+s/t+s^x^t+s/t  ^t+s 


(1.36) 


Finally,  Eqs.  (1.30)  through  (1.35)  feature  the  extended  Kalman 
filter  and  they  can  readily  be  reduced  to  the  Kalman-Bucy  filter  when  the 
systems  are  linear.  Unlike  the  linear  system,  the  covariance  matrices  cannot 
be  precomputed.  And,  as  a  matter  of  fact,  they  are  coupled  with  the  optimal 


estimate  through  coefficients  f  and,  hence,  they  are  not  ordinary  covar¬ 


iance  matrices,  but  rather  they  represent  conditional. covariance  matrices. 


1.4  The  Problem  to  be  Studied 

The  problem  treated  in  the  subsequent  study  is  that  of  estimating 
the  state  of  a  continuous  nonlinear  dynamical  system  (1.10),  influenced  by 
Brownian  motion,  using  discrete  nonlinear  observations  (1.12)  corrupted  by 
an  independent  Gaussian  noise  sequence.  In  the  previous  section,  the  nonlin¬ 
ear  system  is  linearized  and  the  Kalman-Bucy  linear  filter  theory  is  applied 
to  the  problem,  of  estimating  the  state  of  the  linearized  system.  This  tech¬ 
nique  is  based  on  the  assumption  that  the  state  deviation  is  small  so  that 
the  second  or  higher  order  terms  in  the  Taylor  series  expansions  can  be  neg¬ 
lected  while  retaining  the  first  order  terms.  Suppose  that  the  system  is 
highly  nonlinear  or  that  the  initial  uncertainty  is  rexatively  large  so  that 
the  square  of  the  state  deviation  as  well  as  the  deviation  itself  is  not 
negligible.  In  this  situation,  the  Kalman-Bucy  linear  filtering  theory  must 
be  abandoned  and  an  effort  must  be  made  to  develop  a  new  tneory  which  hopefully 
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applies  to  both  linear  and  nonlinear  systems. 

Up  to  the  present  time,  no  sign  of  an  exact  solution  to  the  nonlin- 
'■  !r-  osLitv.tio::  pro!  ,U:::  is.  seen ,  unless  cnc  calculates  an  infinite  number  of 
moment::.  Therefore,  some  sort  of  an  approximate  solution  ii>  ‘‘nevitable.  With 

the  possible  exception  of  a  scalar  system,  it  is  not  practically  feasible  to 

£ 

include  terms  of  higher  order  than  the  second  order  and,  hence,  it  is  desir¬ 
able  that  the  nonlinear  estimation  technique  be  defined  by  using  only  the 
first  two  moments,  namely,  the  conditional  mean  and  the  conditional  covari¬ 
ance.  In  order  to  do  so,  the  second  order  terms  are  included  in  Taylor  series 
expansion  and  a  minimum  variance  criteria  is  employed  to  find  the  conditional 
expectation.  By  definition,  the  conditional  covariance  matrix  is  nothing  more 
than  a  conditional  mean  of  the  square  of  the  actual  estimation  errors  and, 
furthermore ,  ,i.l  is  clearly  understood  that  the  square  of  the  actual  estimation 
errors  js  a  random  Variable.  Therefore,  it  is  meaningful  to  interpret  the 
conditional  covariance  matrix  as  the  optimal  estimate  of  the  square  of  the 
actual  estimation  error  and  to  approximate  it  by  the  same  technique  as  the 
conditional  mean  is  approximated.  It  is  necessary  to  define  a  meaningful 
criteria  for  approximating  the  covariance  matrix  and  this  is  accomplished  by 
the  use  of  the  property  of  the  trace  of  the  matrix.  This  property  is  discussed 
in  Appendix  E. 

The  nonlinear  estimation  theory  developed  in  this  study  is  applied 
to  an  orbit  determination  problem.  The  actual  model  employed  involves  the 
investigation  of  the  states  of  an  interplanetary  space  vehicle  during  the 
planetary  fly-by  and  planetary  orbiter  phases  of  the  mission.  In  the  simu¬ 


lated  study,  Jupiter  is  chosen  as  the  main  body  with  the  Sun  as  the  perturb- 
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1.5  Outline  of  Study 

In  Chapter  2,  the  nonlinear  mathematical  model  which  will  be  studied 
is  discussed  briefly.  On  the  basis  of  a  Martingale  property  which  is  presented 
in  Appendix  D,  an  approximate  nonlinear  estimation  algorithm  is  developed. 

In  the  process  of  developing  the  algorithm,  basic  smoothing  properties  des¬ 
cribed  in  Appendix  B  are  extensively  used  to  manipulate  the  lengthy  algebraic 
relations  and  to  simplify  the  resulting  expressions .  First ,  a  sequential  non¬ 
linear  estimate  is  obtained  and  a  formal  limiting  process  is  used  to  obtain  a 
continuous  nonlinear  estimation  algorithm.  In  the  limiting  process,  the  con¬ 
cept  of  white  noise  as  a  time  derivative  of  Brownian  motion  is  essential.  The 
Brownian  motion  is  treated  separately  in  Appendix  C. 

Chapter  3  is  concerned  with  the  physical  problem  to  be  studied  using 
the  nonlinear  estimation  algorithm  developed  in  Chapter  2.  The  problem  is 
that  of  estimating  the  state  of  an  interplanetary  space  vehicle  during  the 
planetary  fly-by  and  planetary  orbiter  phase  of  a  Jupiter  mission.  The  equa¬ 
tions  of  motion  for  the  spacecraft  are  discussed  briefly  and  expressed  as  a 
set  of  nonlinear  state  dynamic  equations.  Four  kinds  of  observations  are  con¬ 
sidered.  They  are  range,  range-rate,  and  sun-planet  and  star-planet  angles 
as  measured  from  the  spacecraft.  Finally,  computer  programs  for  the  nonlin¬ 
ear  and  the  extended  Kalman  filter  equations  are  described. 

In  Chapter  4,  the  results  of  the  numerical  simulations  are  discussed. 
Several  nonlinear  estimation  algorithms  are  obtained  from  the  modified  Gaus¬ 
sian  second  order  filter  which  is  developed  in  Chapter  2 ,  and  the  modified 
truncated  second  order  filter,  by  neglecting  the  second  erder  terms  in  var¬ 
ious  combinations.  Each  nonlinear  filter  in  conjunction  with  the  extended 
Kalman  filter  is  simulated  with  the  problem  discussed  in  Chapter  3  to 
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determine  the  effects  of  the  second  order  terms,  i.e.,  the  dynamic  second 
order  term,  the  observation  second  order  term,  and  the  Kalman  gain  compensa¬ 
tion  term.  The  Kalman  gain  compensated  filter  obtained-  from  the  modified 
Gaussian  second  order  filter  by  neglecting  the  dynamic  and  the  observation 
second  order  terms  while  retaining  the  Kalman  gain  compensation  term  is  shown 
to  bo  the  best  filter  on  the  basis  of  the  simulations.  The  Kalman  gain  com¬ 
pensated  filter  is  further  examined  through  numerous  simulations . 

A  summary  of  results  and  a  list  of  possible  extensions  to  this  work 


are  presented  in  Chapter  5. 


CHAPTER  2 

A  NONLINEAR  ESTIMATION  ALGORITHM 

2.1  Introduction 

The  state  of  the  dynamic  system  is  assumed  to  evolve  as  the  solu¬ 
tion  of  a  nonlinear  stochastic  differential  equation. 


dxt  =  f(xt,t)dt  +  dBt  ,  t  >  tQ 


(2.1) 


which-  is  expressed  formally  as 

dx,. 


sr  =  «vsi  * ", 


(2.2) 


In  the  above  expression,  f(xt,t)  is  a  n-vector  and  (ut,t  >  tQ}  is  an  n- 
vector,  zero-mean,  white  Gaussian  noise  process  with 


E(utuTf]  =  Qt5(t  -  r) 


(2.3) 


where  Qt  is  an  n*n  positive  definite  matrix  for  any  t  .  Suppose  that  ob¬ 
servations  on  the  state  are  taken  at  discrete  instants  of  time  and  s  mea¬ 
sures  the  time  interval  between  a  certain  point  in  time,  say  t  +  s  ,  and  t 
at  which  the  last  observation  was  made.  Therefore,  s  will  vary  from  zero 
to  the  maximum  time  span  between  two  consecutive  observations.  This  approach 
is  necessary  when  observations  are  not  taken  regularly. 

Let  the  observations  of  the  state  be  of  the  form 


y^  =  h(x.)  t  v.  ,  i  =  1,2,. 


(2.4) 


where  y^  and  h  are  m-vectors ,  and  where  fv^,i  =  1,2,...}  is  an  m-vec- 

tor,  zero-mean,  Gaussian  noise  sequence  with 
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The  covariance  matrix  R  is  an  mxm  positive  definite  matrix  for  any  k  .  Jf 
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It  is  assumed  that  ,a  ,t  >  t  }  and  (v.  }  are  statistically  independent.  :  '■§ 

VO  K  j 

An  extension  to  the  continuous  observation  case  can  be  made  by  simply  replac-  j§ 
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f  " 

ing  v.  with  a  white  Gaussian  noise  v  .  In  this  case  the  function  R,  A 

1  .  .  :  3 

will  have  an  infinite  magnitude.  Since  white  noise  is  formally  modeled  as  ?| 

1 

the  time  derivative  of  Brownian  motion  (see  Appendix  C),  it  is  natural  ,  J| 

to  relate  v.  with  a  white  noise  v  as  follows  'll 

1  t  •  M 

B.  -  3.  ■  % 

v.  =  1+5  1  (2.6)  .  ;| 

i  s  W 

1 

lim  ®t+s  ”  11 

v  =  — - 1  =  ^  (2.7)  | 

t  S-K5  S  ds 

s=t  ;  ^ 

With  these  definitions,  would  be  of  the  form  :  '  g 

1 

**  =  r  (2-8)  |  J 

HI 

I*-:/. 

which  approaches  K  6(s)  as  s  goes  zero.  Denoting  -  'jg 

«•'  "  j 

Y.  =  (y.  ,i  =  1,2, ...t}  (2.9)  j-  .  .•? 

1  t  i  ••& 

t  t  -ch 

I 
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j 

3  that  of  estimating  the  state  x  of  the  dynamical  system  (2.2)  on  the  basis.,  t  'M 

I  t  ;  -Am 

i  .  !  •  § 

S  of  Yt  .  In  particular,  the  desired  estimate  is  the  minimum  variance  esti-  >  ^ 

;  fj 

1  mate  and  the  solution  is  well  known  to  be  the  conditional  expectation  ,  -Jj 

i  !  , 
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I  hlx^/i^J  (5;.  The  details  are  discussed  m  Appendix  l. 
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models  are  nonlinear.  Therefore,  an  approximate  solution  is  inevitable.  It 
is  common  practice  to  linearize  the  system  dynamics  f(*)  and  the  observa¬ 
tion-state  relationships  h(*)  ,  aoout  a  specified  reference  trajectory  and 
then  to  apply  the  Kalman-Bucy  filtering  theory  to  the  linearized  system.  In 
this  chapter,  an  approximate  nonlirear  filter,  which  is  a  modified  Gaussian 
second  order  filter,  is  derived  by  utilizing  Martin^  ;le  properties  (Appendix 
D)  and  a  Taylor  series  expansion  of  f(*)  and  h(*;  about  the  current  op¬ 
timal  estimate,  retaining  the  second  order  terms  in  each  expansion. 

Regarding  the  square  of  the  actual  estimation  errors  as  a  collec¬ 
tion  of  random  variables ,  the  conditional  covariance  is  obtained  by  minimiz¬ 
ing  the  following  risk  function  (see  Appendix  E) 

R(W«)  *  ’  W.)T'  (2-10) 

where  V  is  an  n*n  matrix  and  the  square  of  the  actual  estimation  errors 
and  ^t+syt+c.  is  the  optimal  estimate  of  V^+s  given  Y^+s  ,  which  is  the 
conditional  expectation  of  Vt+g  given  Y^+g  . 

2.2  Apriori  Estimate  ^+s/t 

Integrating  the  state  dynamic  equation  (2.2)  from  t  to  t  +  s  , 


the  state  at  t  +  s  can  be  formally  expressed  as  follows 
t+s  t+s 


x^ ,  =  x_  + 

t+s  t 


f(xT,i)dT  + 


j  u  dt 

it  T 


t+s 


=  +  (*L.  -  + 


t/t  '”t 
t+s 


W  +  Jt  “• 


dt 


(2.11) 


♦  y*t/t.o(*T  ••  k/x)  ■  vdT  (s.i2) 
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The  approximate  expression  (2.12)  is  obtained  by  utilizing  a  Taylor  series 
expansion  of  f(*)  and  truncating  at  the  second  order  terms.  In  this  expres¬ 
sion  t  is  merely  a  parameter  and  s  is  a  variable.  For  the  convenience 
of  notation,  r  in  f(*)  ,  f. (•)  and  f  ^(*)  is  neglected  unless  other- 
wise  staled.  Knowing  that  ,x  >  t  is  independent  of  '{ t  ,  and  (x^  - 
and  (x,t  -  have  zero  conditional  means,  the  conditional  expectation 

^t+s/t  wh*ch  -*r*  the  desired  apriori  estimate  of  x.+s  »  is  obtained  after 
taking  the  conditional  expectation  of  both  sides  of  Eq.  (2.12). 


t+s/t 


t+s 

*t/Wt  <«*,/,>  *  I  f: 


‘  PT/tidl  (2'13) 


Note  that  f(£  .)  ,  f  (£  ..)  and  f  (&  ,  )  are  Y  -measurable  and  the 
i/t  x  x/t  xx  x/t  t 

smoothing  property  3  (Appendix  Bl  can  be  applied.  The  differential  equation 
for  xt+„/.  >  as  a  function  of  s  ,  is  readily  obtained  by  differentiating 
Eq.  (2.13)  with  respect  to  s  .  Since  the  upper  limit  of  integration  is  a 
function  of  s  ,  Leibnitz's  rule  is  applied  and  the  result  is 
d2. 


t+S/t  _  -/£  \  j.  i  jc  (■&  >  .  p 

ds  t+s/t ’  2  *xx  t+s/t'  ‘  t+s/t 


(2.14) 


The  above  differential  equation  is  different  from  that  of  the  extended  KaLr-an 
filter  through  the  inclusion  of  the  dynamic  second  order  term  f  (•)  :  P/2 

AA 

ana  must  be  integrated  in  conjunction  with  ^T+^/t  >  from  t  to  the  instant 
of  a  new  observation,  vising  as  the  initial  condition. 

2.3  Apriori  Conditional  Covariance  Matrix  P  /#, 

T  -f  ?  /  t 

Subtracting  out  Eq.  (2.13)  from  Eq.  (2.12),  the  apriori  estimation 
error:  x  ,  at  t  +  s  ,  is  obtained  as  follows 

l  'S/  v 


Sf  , 


7  :;  •  -  ■  -V;~- — 


*t+s/t  "  xtts  Xt+s/t 


=  (xt-*t/t)f 


.t+s 

u  d 

n  T 


(2.15) 


T+S 

\  ». 

Jt  X 


(8t/t)(iit  - 


<pT  -  >w>*- 


Differentiation  of  Fq.  (2.15)  with  respect  to  s  yields  the  following  dif¬ 
ferential  equation  for  apriori  est  im  at  ion  •  error 

•  %+s7+  -+  i  A  *,/'  v-  1 

ds  -  J  Jf  -  fx^.t+s/.t)xti|/.t  *  2  vxx^fs/t:^  ‘!  -;%#§=  (2,1^ 


%  definition-- 


P  =  V  # 

t+s  t+s/t  t 


'C2V1# 


and  differentiating  Eq.  (2.17)  with  respect  to  s  ,  the  following  relation 


is  obtained 


^*t+s/t 


t.+s/t  ,T  - 

ds  t.+s/t  t+s/t  ds 


(2.18) 


Substitution  of  Eq.  (2,16)  into  Eq,  (2.18)  yields 
dP 

"dF1  =  fx(*t+.s/t)Pt+s  +  |{fxx(^tt3/t);:(PtTs"Pt+s/t)}xt+s/t  +  yL/t 

+  Pf  +s  ^x f 3’t  +s/ 1 1  +  ^t+s/t'^xx^t+s/t^pt+s_Pt.+s/±^  +  ^t+s/t^+s 

Since  t  is  merely  a  fixed  parameter,  .the  conditional  expectation  given  Y 
and  time  derivative  can  be  interchanged.  Therefore,  the . following  expression 
is  obtained  after  taking  the  conditional  expectation  of  both  sides  of  Eq. 
(°.19)  and  incerchanging  the  conditional  expectation  and  the  time  derivacive. 


(2.19) 


‘1 
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*  ywtV*  +  c‘u«.<«/e/rt> 

+  Vs/t^W’  +'  ^Wws'V 


(2.20) 


In  the  above,  the  symmetry  of  the  probability  density  function  is  assumed  and, 

hence,  the  third  order  moment  is  taken  co  be  zero.  The  remaining  terms 
T 

E  fut+s*t+s/t/Yt  ^  can  comPu'te<^  as  follows 


;?%Si’L/t/yt]  *  +  L  E["t+SUT/Yt'dT 


(2.21) 


=  ?  Qt+S 


(2.22) 


Hie  factor  —  comes  from  the  property  of  the  delta  function  (6).  By  the 


same  token 


E'-Xt*s/tut+s/,Yt-'  ~  2  ^t+s 


and,  hence,  Eq.  (2.20)  becomes 


+  VS/tfX«tWt>  +  ?t+S 


(2.23) 


which  is  the  desired  matrix  differential  equation  for  the  apriori  conditional 

A 

covariance  matrix  .  The  initial  condition  is  given  as  a  posterior 

'  A 

conditional  covariance  matrix  Vt/t  * 

2<4  Predicted  Observation 

It  is  apparent  from  Eq.  (2.*’)  ’  .’■he  actual  observation  at  t  4  s 

would  be 


y.^  =  h(xJ.J  )  +  v. , 

■^t+s  t-*5  t+S 

Expanding  h(*)  in  Taylor  series  about  the  apriori  estimate  ^t+s/t  and 
neglecting  the  third  or  higher  order  terms,  the  fo3.1owing  approximate  expres¬ 
sion  for  y  is  obtained 

vTS 

vt«  =  (2-J5> 

After  taking  the  conditional  expectation  of  both  sides  of  Eq.  (2.25)  given 
Yt  ,  the  predicted  observation  ds  obtained. 


^t+s/t  "  h(*t+s/t^  +  2  hxx^t+s/t^  !  Pt+s/t 


(2.26) 


The  above  relation  is  different  from  that  of  the  extended  Kalman, filter  through 
the  second  order  term  in  the  observation-state  relation,  i.e.,  h  (•)  :  P/2  . 

XX 

The  expected  errors  between  actual  and  predicted  observations  are 
obtained  by  the  difference  in  Eq.  (2.25)  and  (2.26).  Hence,  the  apriori  ob¬ 
servation  error  (or  residual)  is  given  by 


t+s  ^t+s/t 


=  h- /*)*. 


+  xh  (*__^)  :  (P_„  -  P. 


x-  t+s/t  t+s/t  2  xx  t+s/t 


)  *  V. 


t+s  t+s/t  t+s 


(2.2?) 


2.5  Posterior  Estimate  and  the  Optimal  Gain 

According  to  the  Theorem  2  of  Appendix  D,  the  following  sequence 

^t+s/1 5^t+s/2  ’*  *  *  ’^t+s/t ’^t+s/t+s ’’  (2*28) 

constitute  a  Martingale  and  if  zi »z2 *z3» •  are  defined  as 


*t+s/i*z2  =  *t+s/2  "  *t+s/l*z3  "  *t+s/3  ~  *t+s/2; 


*t+s/t  "  *t+s/t~l » ; *  ,z^  s  =  *t+s/t+s  '  *t+s/t * "  * 


(2.29) 
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satisfy  the  following  conditions 
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^tlzj  ]  <  “  ,  ^Zn+l/VV‘*,,zn]  =  0  »  n^O 


(2.30) 


with  probability  i.  Further,  the  terms  *t+s/n  ,  n  =  i,2,...  of  the  se¬ 


quence  (2i28)  are  partial  sums  of  the  series  T  z  .  In  other  words,  if 


zt+s  is  so  determined  that  the  conditions  (2.30)  are  satisfied,  the  term 


^t+s/t+s  -°^  se3uenPe  (2.26),  which  is  the  posterior  estimate,  is  uniquely 


t+s 


determined  as  a  partial  sum  V  z  .  Since  z  _  must  be,  a  function  of  6b- 

n=l  x‘rs  ’  -  " 


serration  yt+a  ,  the  following  linear  approximation 


Zt+s  ^t-*s^ t+s  J‘  bt+s 


(2.31^) 


is  assumed.  Where  Kt+s  and  bt+g  are  random  variables  which  are  measur¬ 


able  over  the  o-field  generated  by  the  observations  Y  .  The  bias  term 

.  '  "  ' '  -  t 

bt+s  ^  given , as 


t+s 


i  _K  -  O 

t+s^tfs/t 


from  the  condition  (2.33)  which -must  hold  for  the  z  's  .  From  the  series 

. .  ■  n 


(2.29)  and  Eq.  (2.31) ,~ the  posterior  estimate  ^+s/£+s  is  expressed,  there¬ 
fore,  as  follows 


t+s/t+s  ^f+s/t  +  Kt+s^t+s  "  ^t+s/t^ 


(2.32) 


Tn  Eq.  (2.32),  the  n*m  matrix  K  can  be  chosen  from  a  familv  of  Y  -mea- 


.  4-+s  - -  -  't 

surable  functions  so  that  the  minimum  variance  or  equivalently  the  -minimum 

of  the  trace  of  E[V  J  is  achieved.  By  definition. 
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(2.33) 


=  trE[Cxt+s  -  *t+s/tt8>(*tts  -  *tts/tts>Tl 

trE[{xt+s  -  £t+s/t  ~  Kt+s(yt+s  ~  ^t+s/t^ 

^xt+s  Rt+s/t  ^t+s^yt+s  ^t+s/t^  ^ 
t 

and  the  optimality  condition  (Appendix  F)  for  the  minimum  of  the  above  rela¬ 
tion  requires  that 

E[{(xt*s  ~  *t+s/t^  “  ^t+s^t+s  ”  ^t+s/t^^t+s  "  ^t+s/t'i  //Yt-*  =  0  ^*34^ 

Substituting  Eq.  (2.27)  into  Eq.  (2.34)  and  using  the  smoothing  property  3  of 
Appendix  B ,  the  following  relation  is  obtained. 


E'(  WWi>  V  W>(WW>  + 

*  W«VW(WW ♦  k^W^'WW^m1  <i-35) 

<vw)0ww + k<W'V‘tw*l,’«*,T,r«l 

s' 

If  the  estimation  errors  are  assumed  to-be  jointly  Gaussian,  and  if  it  is 
assumed  further  that  -  - "" 


E«hxx(*t«/t>  !  Pt+=H!,xx(!W<:>  :  Ptis!  '\l 


°  3fhxx(l!t+s/t)  ;  Pt+s/t,{hxx(2t+s/t)  :  W1 


then,  the  optimal  gain  K^+s  is  given  by 


(2.36) 


Kt+s  =  Et+s/thx(ilt+s/t)[hx(at+s/t)^t+s/thx(5it+s/t)  +  Rt+S 
*  i  ^xx^t+s/t5  :  Rt+s/tHhxx(Rt+s/t)  :  W}  ]_1 


(2.37) 


26: 


The  approximation  given  in  (2.36)  leads  to  the  identical  results  obtained  by 
Athans  et  al.  (27)  for  the  scalar  case  but  yields  a  slightly  different  result 
for  the  vector  case.  However,  in  contrast  to  the  result  obtained  by  Athans 
et  al.  the  Kalman  gain  compensation  term 

2*  *hxx^t+s/t^  :  ^t+s/t^hxx^t+s/t^  '  ^t+s/t* 

-  _  -  -  .  —  T  -  ~ 

is  always,  positive  definite,  as  are  the  .terms,  ^xPt+s/t^x  an<^  \+s  °* 

(2i37).  Hence,  the  matrix  to  be  inverted  in  Eq.  (2.37)  will  always  be  posi¬ 
tive  definite,  for  non-zero  .  In  the  extended  Kalman  filter,  the  op¬ 

timal  gain  K  can  be  expressed  in  terms  of  either  apriori  or  posterior 
conditional  covariance  matrix,  which  is  given  in  Eqs.  (1.35)  and  (1.36).  How¬ 
ever,  this  situation  is  not  possible  in  Eq.  (2.37),  due  to  the  Kalman  gain 
compensation  term. 


2.6  Apriori  Estimate  Vt+g/t  of  Vt+g 

By  virtue  of  random  variable  x  ,  the  posterior  estimation  error 

TtS 

*t+s/t+s  *s  an  n><1  vector  random  variable,  and,  hence,  Vt+g  represents  an 
nxn  matrix  of  random  variables  and  its  apriori  estimate  which  is  the  condi¬ 
tional  expectation  of  given  can  be  obtained  from  the  definition 


t+s 


=  * 


t+s/t+s  t+s/t+s 


^t+S  ~  ^t+s/t  ^t+s^t+s  ^t+s/t^ 

Substitution  of  Eq.  (2.27)  into  Eq.  (2.38)  yields  the  following  development, 


Vt+s  '  ‘(I  -  VsVW  -  Wf  hxx  ‘  (Pt«  -  W  +  Vs» 

[(I  -  ’Wx’W  -  Kt«(|  hXX  ’  (Pf.  -  Pt+S/t>  +  Vt«»T 


or 


t+s 


[I  -  ‘W'x'W1  -  Kt«hx1 


.  ,i 


-  <(I  -  Wx'W.1!"*  '  (?t«  -  W>  * 


(2.39) 


Kt+s{2  hxx  :  <Pt+s  “  Pt+s/t5  +  vt+s}*t+s/t{1  “  '^t+shx}T 

where  the  argument  of  h  and  h  is  $  ,  . 

X  XX  T+S/T 

Knowing  that  v^+g  is  independent  of  Yt  and  taking  the  condi¬ 
tional  expectation  of  both  sides  of  Eq.  (2,39)  given  Yt  yields  the  follow¬ 
ing  approximate  expression  for  the  apriori  estimate  of  V  , 

L  »  5 

ft+s/t  '  {I  '  Kt«ViW>>W<I  -  Kt«h*<*t«ft»T 

(2.40) 

+  Kt+8^7  ^XX^t+s/t^  :  Pt+s/t){hXX(\+s/t)  :  Pt+S/t}  +  ^t+S^t+S 


?*7  posterior  Conditional  Covariance  Matrix  V 


t+s/t+s 
>n, 

random  variables.  The  conditional  expectation  V 


As  pointed  out  in  the  previous  section,  V  is  a  collection  of 

t+s 


t+s/t+s  glVen  Yt+S  1S  a 


posterior  estiniate  of  V^+g  and  can  be  determined  by  estimating  each  element 

of  vtis  terns  of  the  linear  combination  of  the  apriori  estimate  V 

-  t+s/t 

and  a  new  observation  yt+g  .  The  above  argument  is  based  on  the  Martingale 
properties  of  the  sequence  \+s/1  .^t+s/2  >•  •  •  >Vs/t , . . . ,  and,  hence,  the 
sa'me  reasoning  as  that  used  in  regard  to  Eq.  (2.31)  can  be  applied.  Regard¬ 
ing  Vt+s  as  an  n?*l  vector  instead  of  n*n  matrix,  the  estimation  problem 


can  be  stated  as  follows:  determine  n^*m  matrix  B 

t+s 

bination 
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of  the  linear  com- 


Vt+s/t+s  "  Yt+s/t  +  Bt+s(yt+s  '  ^t+s/t5 
such  that  the  risk  (see  Appendix  E) 

R(B.  t  )  =  trE [ ( V  -  V  W  xT, 

t+s  t+s  t+s/t+s^'Vt+s  t+s/t+s^ 


(2.41) 


(2.42) 


is.  minimized.  Substituting  Eq.  (2.41)  into  Eq.  (2.42)  leads  to  the  following 


result 


R(B.  ,)■  =  trE[{V. ,  -  V,  .  -  b  (v  -  v  M 

vt+s  —  --  .t+s  -  t+s/.t  t+s-  yt+s  yt+s/t"* 

sm  a  2  *n  (2;43) 

t+s  t+s/t  “  Bt+s^yt+s  “  yt+s/t)}  J 

The  optimality  condition  for  the  minimum  of  R  is  given  by  the  following  or¬ 
thogonality  condition  (see  Appendix  F) 

E[{Vt+s  ^t+s/'t  Bt+S^yt+S  "  ^t+s/t^yt+s  ~  Wt^V  =  0  (2.44) 

E  [(vt+s  "Vt+s/t }  (yt  ts^t+s/t )T/ Yt ]  =  Bt+SE  [  <yt  +sA+s/t }  (yt+s"^t+s/t  )/Yt3  (2‘45) 

By  the  same  argument  as  that  used  with  regard  to  Eqs.-  (2.34)-,  (2.35)  and 
(2.36)  the  right  hand  side  (R.H.S.)  of  Eq.  (2.45)  is  approximated  as  follows 

t+sl  t+s/t'  t+s/tV*t+S/t'  +  Rt+S 


+  7  (hK*(W>  >  <  ft+s/t>Tl  <2'46> 


Knowing  that 


E^t+s/t^yt+s  "  ^t+s/t^  ^Yt^  =  0 


(2,47) 


the  left  hand  side  (L.H.S.)  of  Eq.  (2.45)  becomes 
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L.H.S.  =  E[Vtts:yt+s  -  ?t+s/t>  /Yt) 


(2.48) 


After  substitution  of  Eqs.  (2.27)  and  (2.39)  into  E^.  (2.48),  it  can  be  easily 
seen  that  all  of  the  terms  in  Eq.  (2.48)  are  of  the  fourth  or  higher  order 
moments,  under  the  assumption  that  the  probability  density  functions  of  es¬ 
timation  eri'ors  are  jointly  Gaussian.  Therefore,  the  optimal  gain  B^+s  is 

given  as  the  ratio  of  Jthe  fourth  to  the  second  order  moments  and  is  neglected. 

/ 

Hence , 


]1  =  0 
t+s 


(2.49) 


With  this  assumption,  Eq.  (2.L1)  l-ads1  to 


(b  =  0 

t+s/t+s  t+s/t 


(2.50) 


It  is  interesting  to  note  that  both  the  posterior  and  the  apriori  estimation 
errors  are  independent  of  observations  available  for  the  linear  model  and, 
hence,  the  conditioning  on  the  covariance  matrices  becomes  unconditional. 
Therefore,  there  are  no  uis+inctions  between  ^t+s/t+s  and  ^ t+s/t  ’ 
though  Eq.  (2.50)  shows  that  Vt+s/t+s  is  closely  approximated  by  Vt+s^t  , 
these  are  conceptually  two  different  quantities.  In  the  linear  filtering 
theory,  these  become  identical  and  there. is  no  distinction  between  them. 

Since  is  related  to  the  apriori  conditional  covariance 

matrix  through  Eq.  (2.40),  the  posterior  conditional  covariance  matrix 

^t+s/t+s  can  k®  expressed  in  terms  of  •  From  Eq,  (2.50) 


V  ~  v 

t+s/t+s  t+s/t 


•  (I  -  1W,x(JWl:))ftts/t{I  - 


t+S  X  t+s/t' 


(2,51) 


+  Kt+s  ^2  *hxx**t+s/t*  :  ^t+s/t^xx^t+s/t5  :  Rt+s/t*  +  Rt+s^Kt+s 


pi 


P:  I 

Plf*>  m 


gp 

hi 

f§ 

V%-  tP 


I  «::/ 1  t:. 


'  t  )■:/.  •  \1„\<‘l.Wl:)W  •  Plt,/t"IW.'«!/1,KL 


+  K  fh  (#  .  .  )?.  ,  .,h  (x.  ..)  t  R 

1:i-:  x  t+s/t  t+s/t  x  t+s/t  t 


(2.52) 


•  2  '  XX  "t+s/t'  *  t+s/tJl  xx^t+s/t'  •  t+s/t J  J1\+s 

Substitution  of  Eq,  (2.36)  into  Eq.  (2.52)  yields  the  following  relationships 


W«  ■  W  -  W’Wt  (2-53) 


The  above  relation  is  used  to  update  the  apriori  conditional  covariance  matrix 

a  j% 

Rt+s/t  t0  posterior  conditional  covariance  matrix  Vt+s/t+s  after  a  new 
observation  yt+s  is  processed.  Once  \+s/-j-+s  Stained,  it  can  be  used 
as  an  initial  condition  for  the  integration  of  Eq.  (2.23)  from  t  +  s  to  the 
instant  of  a  new  observation.  Finally,  the  procedures  required  to  compute 
the  posterior  estimate  ^ts/tts  can  summarized  as  follows: 

=  «w>  ♦  \  yw.1 :  w  <2-54) 


r 

^-ts/t  .  £  /A  \p  .  p  )  +  0 

ds  3i  t+s/t  t+s/t  t+s/t  x'  t+s/t  ^  Vs 


(2.54) 

(2.55) 


^ t+s/t  =  h(V../J  :  K 


t+s/t  2  xx  t+s/t  *  t+s/t 


(2.56) 


Kt+S  =  ^t+s/thxUt+s/t)Ihx(^t+s/t)^t+s/thx^t+s/t)  +  Rt 


+  £<***<*++./+> ? 


(2.57) 


r  2  l"xx'  t+s/t  *  1  t+s/t xx^t+s/t'  ’  1  t+s/t J 
t+s/t+s  ~  ^t+s/t  +  Rt+s^t+s  “Vs/t^ 


V  s  p  -  K  h  (ft 

t+3/t+S  t+s/t  t+SXV  t+s/t  t+s/t 


(2.58) 

(2.59) 


In  order  to  start  the  computation,  %Q/0  a^d  Vq^o  are  required  and  they 
are  given  by 


I 


3i 


x  ,  =  E[x  n  3  =  E[x  ] 

o/o  1  o  o  1  o 


(2.60) 


V  .  =  P 

o/o  o 


(2.61) 


where  P  is  given  as  the  covariance  cf  the  random  variable  x  . 
o  o 

Depending  on  the  particular  problem,  it  may  be  possible  to  neglect 
either  the  dynamic  or  the  observation  second  order  term.  This  might  be  the 
case  when  the  state  dynamics  are  relatively  smooth  while  the  state-observa¬ 
tion  relationships  are  highly  nonlinear  or  vice  versa. 


2.8  Computational  Algorithm 

The  algorithm  for  computing  the  estimate  9  .  of  x  by 

t+s/t+s  t+s 

processing  each  data  point  sequentially,  can  be  summarized  as  follows: 

1.  Given  «t/t  and  «t/( 

2.  Compute  *t+s/t  an*3  ^t+s/t  integrating  Eq.  (2.54)  and 

(2.55)  with  the  given  initial  conditions  9  and  V  . 

t/t  t/t 

until  a  new  observation  v^,  is  made  at  t  +  s  . 

"  t+s 

-  3.  Determine  ^t+s/t  and  Kt+s  using  Eqs.  (2.56)  and  (2.57), 

respectively. 

4.  Compute  *t+s/t+s  and  $t+s/t+s  by  updating  *t+s/t  and 
Pj+s/t  through  Eq.  (2.58)  and  (2.59),  respectively. 

5.  Given  ^t+s/t+s  and  ^t+s/t+s  >  the  steps  2  through  4  car.  be 
repeated. 


2.9  Continuous  Second  Order  Filter 

An  approximate  filter  for  the  case  of  continuous  observation  may 
be  obtained  by  passing  to  a  formal  limit.  In  doing  so,  Qt+S  the  dynamic 
state  noise  and  R^+g  ° f  the  observation  noise  have  to  be  replaced  by  Q^s 


and  R^/s  ,  respectively.  This  comes  from  the  property  of  white  noise  re¬ 
garded  as  a  time  derivative  of  Brownian  motion  (see  Appendix  C).  For  an  in¬ 
finitesimal  observation  interval  s  ,  it  follows  that 


t+s/t 


2  +  sf(*,J  :  V 


't/t 


t/t'  T  2  *xx^t/t'  ‘  vt/t 


(2.62) 


U  5  Vt +  sfxwt/t^t/t +  ^t/tfXt)  +  sQt  (2-63) 


K 


-  T/  -  \1J-1 


t+S 


sPtWtV’Wt)V 


(2.64) 


Utilizing  Eqs.  (2.56),  (2.58),  (2.59),  (2.62),  (2.63)  and  (2.64)  and  passing 
to  a  formal  limit,  the  following  continuous  second  order  filter  is  obtained, 
as  s  goes  to  zero. 


dft 

dt 


o  -  ft 

lim  t+s/t +s  t/t 


s+o 


dft 

dt 


f (*>  +  ~  fxx(ft)  s  V  +  K[y  -  h(ft)  -  -  hxx(ft)  :  V] 


K 


^hx(x)R~1 


(2.65) 

(2.66) 


dV 

dt 


lim  \+s/t+s  Vt/t 


s->-o 


—•  -  f  (ft)V  +  VfT(ft)  +  Q  -  VhldOK- '-hVOV 

at  XX  xx 


(2.67) 


Note  that  the  optimal  gain  K  for  the  continuous  filter  is  not  the  limit  of 


the  optimal  gain  K 


t+s 


given  in  Eq.  (2.64). 


CHAPTER  3 


DESCRIPTION  OF  THE  ORBIT  DETERMINATION  PROBLEM 

3.1  Introduction 

In  order  to  compare  the  performance  of  the  extended  Kalman  filter 
and  the  various  forms  of  the  nonlinear  second  order  filters  developed  in 
Chapter  2,  the  methods  are  compared  in  a  simulated  study  of  a  realistic  or¬ 
bit  determination  problem.  The  problem  considered  is  that  of  estimating  the 
state  of  an  interplanetary  space  vehicle  during  the  orbiting  and  planetary 
fly-by  stages  of  a  Jupiter  exploration  mission.  The  reason  for  choosing  this 
problem  is  that  considerable  attention  has  been  given  to  the  exploitation  of 
deep  space  and  the  reconnaissance  of  Jupiter  is  regarded  as  an  important 
scientific  objective.  However,  the  past  Jupiter  encounter  missions  are  of 
comparable  significance  and  those  missions  are  made  practical  by  utilizing 
the  powerful  trajectory  shaping  capabilities  of  Jupiter's  gravitational  field. 
One  such  mission,  the  so-called  "Grand  Tour",  involves  successive  flv-bys 
of  the  planets,  Jupiter,  Saturn,  Uranus,  and  Neptune.  The  Grand  Tour  is  the 
subject  of  considerable  current  interest,  since  a  mission  opportunity  occurs 
in  the  last  half  of  the  IBVO's  and  will  not  reoccur  for  another  179  years. 

A  critical  problem  in  ti.n  ^«sign  of  a  space  vehicle  to  perform  a 
deep  space  mission  such  as  the  Grand  Tour  is  the  accurate  determination  of 
the  expected  trajectory  which  is  the  basic  knowledge  required  for  the  guidance 
correction.  Because  of  numerous  sources  of  error,  the  true  trajectory  is 
never  known  to  us.  A  major  contribution  to  those  errors,  access  during  the 
encounter  trajectory,  due  to  imperfect  pre-encounter  guidance  corrections 
which  rosult  from  pre-encounter  orbit  determination  errors.  In  regard  to 
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this  situation,  it  is  interesting  to  see  how  the  nonlinear  estimation  proce¬ 
dures  perform  when  compared  with  the  extended  Kalman  filter. 

3.2  Equations  of  Motion 

The  motion  of  a  space  probe  relative  to  a  given  planer  is  closely 
approximated  by  the  solution  of  the  following  vector  stochastic  differential 
equations 


-§r 


-  1'  l  jr  -  -rl  +  U 


(3.1) 


where  y  and  y  are  the  gravitational  constants  of  the  target  planet  and 
s 

the  Sun,  respectively,  and  u  is  a  vector  of  Gaussian  process  noise  and 
where  r  is  the  position  vector  of  the  space  probe  relative  to  the  target 
planet,  r^  is  the  position  of  the  target  planet  relative  to  the  Sun,  and 
r^  =  r  +  is  the  position  of  the  probe  relative  to  the  Sun. 

Eq.  (3.1)  can  be  reduced  to  a  system  of  first  order  differential 
equations  by  the  following  transformation 


r  =  v 


(3.2) 


-  ps  (pr  -  jrl  *  « 
P  t 


In  a  cartesian  coordinate  system  centered  at  the  target  planet,  the  equations 
of  motion  can  be  expressed  in  component  forms  as 

X  *  U 


0  m 
p  & 


Y  -  V 


Z  =  W 


(3.3) 


tv  X  +  X 

•  X-  .  t  t, 

u  :s  : -v  ^  -  Ms  [~pr-  -  pT:  +  ux 

•/  p  t 

4.  Y  Xt+Y  Yt 

v  =  -ril  p-  -  ug  .[— --3- -■  -  —yj  +  U 

-  ■  .-  •  P  t 

.,: ;  z.  rzt  +  z  zt 

V‘  .=•  ~u  .£y  ~  I — r3 - rrl  +  uz 

P  t 


(3.3) 


few**  •  '  . r 

y .  •  A.?1. 


Th'e  heliocentric  posit ion  components  of  the  target  planet  can  be  expressed  as 


X-  =  X*  +  b 
<t  .  t  x 


Y:.  -  Y*  +  b 

t  t  y 

Z.  =  Z*  +  b 

~-X  t  z 


PH  ■■  -£  :  :r  :  -  .•:  ,z.  =  z*  +  b 

I  -W-v,  >;.•••  •  .:•  .  .  .r.  •;  t  t  Z 

VS  ^  14.  X*  ,,  Y*  and  ,Z*  are  the  components  of  the  heliocentric  position  vec- 

f-:-V  AS'-  L  *.  ''  J'  ' 

!  'V;  ^brs;;;pbtaxned  fr.of^tbe  planetary  ephemeris  and  b  -  ,  b  and  b .  are  com- 

?.  -2-  V-'v--  ■■,  Vj  -  'X  y  Z. 

.\  «•  sag'  ; 

I  ?  ponents  of^bias  inythe  position  vector  due  to  the  errors  in  the  planetary 

:  ‘  '  !%?w  - 

£  ,  y  ephemeris . '  The  cqmpphents  of  the  planetary  position  bias  are  assumed  to  be 

}  ^  overthe  timeyperiod  of  interest.  The  position  vector  of  target 

'pi#11^  changes  very  siowly  and  this  assumption  appears  to  be  reasonable. 


|v" 


'3.3  Augmented  State  Vector 

Since  the  uncertainty  ir.  the  position  of  the  outer  planets  (Jupiter 
ip  this  study)  is  assumed  to  be  an  influential  error  source,  the  planetary 


bias  vector  b^  ,  b^.  and  bz  are  assumed  to  be  unknown  parameters  and  are 

e?ti?nPPi^-'  fo  achieve  this  objective,  the  original  state  vector  given  by  Eq. 

(-3..  3)  expanded  to  include  b„  ,  b  and  b  and  the  augmented  state  is 

.  y  z 

.  defined  as  -  _ 
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XT  =  [•r.JY,z,u.,v,w,bx,by,bz] 

The  augmented  state  vector  is  governed  by  the  ( 9x1 )  vector  differential  equa¬ 
tion 

(3.5) 


(3.6) 


3.4  St  ate -Ob  serve  t ion  Relationships 

There  are  four  types  of  observations  considered:  range  (p)  , 
range-rate,  -(p)  ,  sun-planet  angle  (a)  ,  and  star-planet  angle  (3)  . 

The  first  two  of  these  are  Earth-based  while  the  other  two  are  onboard  obser- 

/ 

vatiohs.  Any  combination  of  the  above  four  observations  can  be  processed 

at  any  time  interval.  Such  a  procedure  is  necessary  if  the  characteristics 
of  each  type  of  observation  are  to  be  determined. 

1.  The  range  measurement  is  given  by 

P  =  (p  •  P)1^2  +  vp  (3.7) 


which  in.  component  form  becomes 

a  =  cos"1  [(XXt  +  YYt  t  ZZ±  +  r2)/rrp)  +  va  (3.13) 

where 

.  "•  r  =  [X2  +  Y2  +  Z2]1/2 

and 

r  =  [X2  +  Y2  +  Z2]1/2 

p  P  P  P  -  - 

_  4.  Star-planet  angle 

The  last  star-planet  angle  measurement  0  ,  defined  as  the  smaller 
angle  between  the  probe-planet  line  and  the  line  from  the  probe  to  a  reference 
;$gar^  .asl  given  by 

v.  ’  „  C 

0  -  =  COS-1  +  Vg  (3.14) 

where:  is 'the  random  observation  error  and  S  is  a  unit  vector  in  the, 

^  •  •  .  p  -y  -  ■  . 

.direction  of  the  reference  stab.  The  star  is  assumed  to  be  at  an  infinite 
distance  so  that  S  is  a  .constant  vector.  Since  the  inclination  of  Earth 

t 

and  Jupiter  €0*6  nearly  zero,  it  follows  that  Earth,  Jupiter  and  the  space 
probe  lie  in  very  nearly  the  same  plane.-  Hence,  it  is  desirable  to  use  a 
-'star  which'  is  not  in  this  plane  as  a  reference  star  to  obtain  information  - 
about  the’  out-of-plane  motion. '  Star-planet  angle  measurement  can  be  ex¬ 
pressed  in  cartesian  components  as 

8  =  cos”1  (-(XS  +  YS  +  ZS  )/r]  +  v.  (3.15)- 

x  y  z  p 

where  S  ,  S  and  S  are  the  direction  cosines  of  the  reference  star 
x  y  z 


direction. 


,3r!i  ^Motion  of  theTrackihg;  Station 

Range,  .and-  range-rate  observations  are  taken  at  a  tracking  station 
on  the  .Earth  and  the  heliocentric  position  vector  of  the  tracking  station 


is  given  by 


r.  =  r  +  R 
s.  e 


(3.16) 


;where  -f-  is'  the  heliocentric  .position  of  the  Earth  and  R  is  the  geocen¬ 
tric  position  of  the  tracking  station.  The  vector  *R  is  computed  :as.  a 
function  of  time  from  the  relationship  ’  • 


Tr<  cos  6  cos  a-  (t) 

I  '  S  S': 


R  =-  [T]  U  cos  5  sin  a 

-  s  -  _s 


I R  sin  6 , 


-  (3.17)' 


where  ag  is  the  right  ascension  of  the  tracking  station, 

6  is  the  declination  (latitude)  of  the  tracking  station, 

/  5  r  -  - 

■'  j  V  R  is  the  magnitude  of  the  vector:  iR<^aridr-the  .radial  distance^of  the 
tracking. station  from  the  Earth’s center. 


0  cos  e  sin  e 


0  -sin  e  cos  e 


is  the  rotational  matrix  which  transforms  the  coordinates  from  an  equatorial 
to  an  ecliptic  cooroinate  system  which  is  chosen  to  be  the  heliocentric  sys¬ 
tem.  The  argument  e  is  the  obliquity  of  the  ecliptic.  With  the  assumption 
that  the  Earth’s  rotation  is  uniform,  the  right  ascension  of  the  tracking 
station  can  be  expressed  as 


Tracking  Station 
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a.Ct)  =  VV  +  *s  •  (t  -  tQ) 

Differentiation -of  Eq.  (3.16)  with  respect  to  the  time  yields  the  heliocen¬ 
tric  velocity  of  the  tracking  station. 


r . 
s 


re  +  t 


R  cos  6  sin  a  (t) 
s  s 

R  cos  6  cos  a  (t) 
s  s 


Lo 


(3.18) 


The  zenith  angle.  <j>  of  the  probe  with  respect  to  the  tracking  station  is 
.  -given  by 


cos  •  4> 


R  •» .  p 

Rp.,r' 


(3.19) 


The  probe  is  assumed  to  be  visible  from  the  tracking  station  if  cos  <p  is 
positive. 

3.6  Simulation  of  Errors  ~ 

Each  component  of  the  noise  in  the  equations  of  motion  (3.1)  and 
observation-state  relations  is  modeled  as  a  normally  distributed  scalar  ran¬ 
dom  variable  with  zero  mean  and  known  variance .  The  noises  are  simulated 
by  sampling  ait  random  from  a  standard  normal  distribution  function  (zero 
,mean  and  unit  variance)  and  then  scaling  the  sampled  number  by  the  given 
standard  deviation. 

The  normal  density  function  of  the  random  variable  K  is  given  by 


ftt)  s  ~~r  exp  [-  --"hr-]  (3.20) 

<2*0  20 

where  m  and  a  are  the  mean  and  the  standard  deviation,  respectively.  Eq. 
(3.20)  can  be  written  in  terms  of  the  standard  normal  distribution  function 


by  the  transformation 


For  the  random  variable  of  zero  mean 


(3.22) 


a? 


(3.23) 


Ihe  inverse  of  Eq.  (3.21)  can  be  approximated  by  the  curve  fit  equation  (40) 

z  2  r  - 


c  +  c  r  +  c„r2 

o  1  2 

i  +  d.r  +  d„rk  +  d^r* 


(3.24) 


r  =  [*n(F"2)]2 


where  the  coefficients 

c.  and  d. 

1  X 

have  the  following  values 

0 

o 

=  2.515517 

yx 

II 

1.432788 

C1 

=  0.802C53 

II 

CM 

rO 

0.189269 

C-3 

=  0.010328 

d3  " 

0.001308 

Sampling  of  the  standard  normal  distribution  is  accomplished  by  first  samp¬ 
ling  at  random  from  a  uniform  distribution  to  obtain  a  value  for  F(0  <  F  <  1) 
and  then  computing  the  standard  normal  random  number  z  from  Eq.  (3.24). 

The  simulated  noise  is  then  computed  as  the  product  of  the  sampled  value  z 

of  the  standard  normal  random  variable  £  and  the  standard  deviation  0  by 
Eq.  (3.23). 

Observational  data  are  simulated  by  adding  random  numbers  £  which 
are  generated  in  the  manner  described  above  to  the  observation  value  computed 
from  the  true  state  and  state-observation  relationships  discussed  in  Section 


+  V 


(3.25) 


Noise  in  the  equation  of  motion  is  simulated  in  the  same  way  as 
above  and  then  added  to  the  equation  of  motion  at  discrete  points  of  time 
which  correspond  to  the  integration  step. 

3.7  Computer  Program  Description 

The  program  NONSTEP  (NONlinear  STate,  Estimation  Program)  is  dev¬ 
eloped  for  comparison  of  the  extended  Kalman  filter  and  the  nonlinear  es- 
timation  algorithms  by  applying  each  to  the  study  of  an  interplanetary  orbit 
determination  problem.  Special  emphasis  is  given  to  the  planetary  fly-by 
mode  although  the  planetary  or  biter  is  considered  also. 

The  program  was  written  in  FORTRAN  IV  for  the  CDC  6600  computer 
system  at  The  University  of  Texas  at  Austin.  Since  this  computer  has  a  single 
precision  word  length  of  sixty  bits,  single  precision  arithmetic  was  consid¬ 
ered  to  be  adequate  for  most  calculations.  The  initial  frame  of  the  program 
was  founded  on  the  existing  program  STEP  (STate  Estimation  Program)  developed 

by  Jones  (28)  at  The  University  of  Texas  at  Austin. 

The  three  basic  functions  of  the  program,  i.e.,  simulation,  estima¬ 
tion,  and  evaluation,  are  conducted  sequentially  according  to  a  schedule 
specified  in  the  input  data.  The  program  NONSTEP  has  a  capability  for  carry¬ 
ing  out  the  nonlinear  estimation  algorithms  as  well  as  the  extended  Kalman 
filter,  depending  on  the  input  data  IFILTER,  If  IFILTER  =  1,  the  extended 

Kalman  filter  is  carried  out.  If  IFILTER  =  2,  the  nonlinear  filter  is 

implemented  and,  finally,  if  IFILTER  =  3,  the  nonlinear  estimation  proce¬ 

dure  is  first  carried  out  and  then,  with  the  same  input  data  and  random 
noises,  the  extended  Kalman  filter  is  carried  out.  With  this  latter  option. 


i*G 


a  direct  comparison  of  the  linear  and  nonlinear  algorithms  can  be  made. 

In  order  to  reduce  the  storage  space  for  compilation,  the  program 
N0N3TEP  employs  OVERLAY.  The  main  OVERLAY(0,0)  controls  the  overall  program. 
The  0VERLAY(1,0)  conducts  all  the  plots  of  the  output  data  for  the  case  of 
IFILTER  =  1  or  2.  The  OVERLAY(2,0)  does  the  same  thing  for  the  case  of 

IFILTER  =  3.  The  0VERLAY(3,0)  conducts  all  of  the  calculations  involved  in 

simulations,  estimations,  and  evaluations  and  transfers  the  output  data  to  a 
magnetic  tape  for  the  plot  of  OVERLAY (1,0)  or  0VERLAY(2,0). 

In  conjunction  with  the  apriori  conditional  covariance  matrix,  the 
true  trajectory  and  apriori  estimate  are  generated  simultaneously  through 
parallel  numerical  integrations  of  the  apriori  estimate  and  the  apriori  con¬ 
ditional  covariance  matrix.  A  general  purpose  numerical  integration  subrou¬ 
tine  is  used  to  simultaneously  integrate  the  differential  equations.  The 
routine  consists  of  a  Fourth  Order  Adams  predictor-corrector  scheme  with  a 
Runge-Kutta  starter.  Although  the  integration  is  carried  out  in  single-pre¬ 
cision,  the  dependent  variables  are  carried  in  double-precision  to  minimize 
round-off  errors. 

Observational  data  are  simulated  by  generating  ra.  lorn  noise  and 
superimposing  it  on  the  true  observation  computed  from  the  true  state, 

A  simplified  block  diagram  of  the  computational  logic  is  shown  in 


Fig.  5. 


CHAPTER  4 


DISCUSSION  OF  NUMERICAL  RESULTS 

The  purpose  of  this  chapter  is  to  determine  the  characteristics  of 
second  order  filters  on  the  basis  of  a  numerically  simulated  study.  There 
are  two  basic  classes  of  second  order  filters  to  be  examined.  The  first  is 
the  modified  Gaussian  second  order  (MGSO)  filter  and  the  other  is  the  modi¬ 
fied  truncated  second  order  (MTSO)  filter.  The  basic  difference  between 
these  two  filters  is  that  the  Kalman  gain  compensation  (KGS)  term  enters  with 
a  plus  sign  of  one-half  in  the  first  filter  and  with  n  minus  sign  of  one- 
fourth  in  the  later  filter.  Both  filters  include  a  dynamic  second  order  (DSO) 
term,  (f  :  P)/2  ,  in  the  dynamic  equation  (2.54)  and  an  observation  second 
order  (0S0)  term,  (h  :  P)/2  ,  in  the  predicted  observation  equation  (2.56). 

4.1  Various  Simplified  Forms  of  Second  Order  Nonlinear  Filters 

Although  the  modified  Gaussian  second  order  filter  and  the  modified 
truncated  second  order  filter  are  developed  using  a  model  in  which  both  state 
and  the  observation  equations  are  nonlinear,  there  is  a  possibility  that  the 
actual  problem  will  consist  of  a  highly  nonlinear  dynamic  equation  and  a  rela¬ 
tively  linear  observation  or  vice  versa.  An  this  situation,  the  second  order 
term  in  the  relatively  linear  relation  may  be  neglected,  and,  hence,  several 
possible  simplified  nonlinear  estimation  algorithms  can  be  obtained,  depend¬ 
ing  on  the  presence  of  the  dynamic  second  order  term,  observation  second 
order  term,  and  the  Kalman  gain  compensation  term  in  various  combinations. 

The  resulting  filters  are  tabulated  in  Table  1.  The  Filters  1,  8  and  9  of 
Table  1  correspond  to  the  MSSO  Filter,  the  MTSO  Filter  and  the  extended 
Kalman  (EK)  Filter,  respectively.  The  Filter  4  is  specifically  referred  to 


as  the  Kalman  Gain  Compensated  (KGC)  Filter.  The  performance  of  each  of  the 
Filters ,  1  through  8 ,  is  studied  through  numerical  simulations  and  compared 
with  the  EK  Filter  which  is  the  most  popular  filter  at  present  time.  The 
nonlinear  filter  is  first  executed  with  the  input  data  given  in  Tables  5 
through  9,  and  then  the  EK  Filter  is  executed  under  the  same  corditions. 

The  same  sequence  of  random  numbers  is  used  to  simulate  the  state  noises  as 
well  as  the  observation  noise  in  both  filters.  The  conclusions  reached  in 
this  investigation  are  based  on  the  results  of  several  hundred  simulations. 
The  results  obtained-  in  eighteen  of  these  simulations  are  .presented  in  this 
report.  These  results  obtained  in  these  cases  are  representative  of  the 
results  obtained  An  the  remaining  studies. 


4.2  The  Nominal  Trajectory  and  Error  Sources 

.  ..  The- nominal  trajectories  are  generated  by  integrating  the  equations 

of  motion  (3.1)  with  dynamic  state  noise  set  equal  to  zero.  In  other  Words,- 
the  random  noise  u  is  simply  set  to  zero.  Considering  the  possibilities 
of  applying  the  nonlinear  algorithm  to  the  problems  of  a  near-Earth  or  lunar 
satellite.  Mariner  and  Viking  missions,  simulations  are  conducted  not  only 
for  a  hypobolic  orbit,  but  also  for  an  elliptic  orbit.  The  nominal,  trajec¬ 
tory  of  the  elliptic  orbit  is  shown  in  Fig.  6  and  the  hyperbolic  orbit  in 
Fig.  7.  The,  periapsis  and  apoapsis  for  the  elliptic  .orbit  occur  at  about 

4 .‘6  and  20  days,  respectively.  The  periapsis  encounter  in  the  hyperbolic 

.  t  -  ,  ’  ' 

orbit  occurs  at  about  12.7  days.  -X' 

___  ■<:  .  - 

As  seen  in  Figs.  6  and  7,  the  dynamic  nonlinearity  in  the  elliptic 

orbit  is  very  much  concentrated  near  periapsis  and  apoapsis.,  but  it  is  well 
distributed  over  the  entire  trajectory  when  compared  with  the  hypobolic  or¬ 
bit.  In  the  hypobolic  orbit ,  the  dynamic  nonlinearity  is  concentrated  almost 
entirely  near  the  perigee,  and  the  pre-  and  post-encounter  trajectories  ap¬ 
pear  to  be  straight  lines. 

The  initial  conditions  for  the  hypobolic  and  elliptic  orbits  are 
obtained  from  the  nominal  Grand  Tour  mission  trajectory  with  Earth  launch 
date  and  Jupiter  encounter  (28).  For  the  elliptic  orbit,  the  velocity  com¬ 
ponents  are  reduced  so  that  it  yields  an  elliptic  orbit  with  a  proper  period 
of  30  days.  The  nominal  trajectory  initial  conditions  are  given  in  Table  2 
and  the  orbital  elements  of  Jupiter  and  Earth  are  listed  in  Table  3. 

The  true  trajectory  (or  state)  is  generated  by  adding  a  vector 
Gaussian  random  forcing  term  u  described  in  Section  3.6,  to  the  equations 
of  motion  (3.1). 


5? 


The  actual  observations  are  simulated  with  the  data  given  in  Table 
4  by  using  the  following  procedure: 

1.  Compute  the  nominal  observation  through  the.  state-observation 
relationship  with  the  true  trajectory  obtained  as  described  above. 

2.  Gaussian  random  noise  is  generated  as  described  in  Section  3.6 
and  added  to  the  nominal  observation. 

The  standard  deviation  Oq  of  the  dynamic  state  noise  a  and 

o_  of  observation  noise  v  are  given  in  Tables  5  through  9  according  to 
K 

the  simulations.  Since  it  is  common  practice  to  employ  an  adequate  o«  for  the 

dynamic  noise  u  ,  even  though  there  is  no  dynamic  noise  assumed,  two  values 

of  an  are  used,  o  is  designated  for  the  true  trajectory  and  on  is 

y  yi  yn 

adopted  for  the  estimate  of  the  state. 
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4.3  Characteristics  of  the  Filters 

The  results  of  Simulation?"  1  through  8  are  presented  primarily  for 
the  purpose  of  describing  the  characteristics  of  the  filter  algorithms  in 
Table  1. 

For  each  simulation,  position  and  velocity  estimation  errors  are 
plotted  and  they  are  compared  with  results  obtained  with  the  extended  Kalman 
(EK)  filter  which  is  filter  9  in  Table  1.  In  addition  to  the  estimation  er¬ 
rors,  the  conditional  variances  and  the  observation  residual  which  is  defined 
as  the  difference  between  the  actual  observation  and  the  predicted  observa¬ 
tion  are  shorn.  Unless  stated  otherwise,  ail  figures  are  obtained  by  con¬ 
necting  every  third  data  point  with  straight  lines.  There  are  thirty  data 
points  between  two  adjacent  symbols.  The  main  reason  for  sampling  every 
third  data  point  is  due  to  the  difficulty  of  tracing  the  original  plot  ob¬ 
tained  from  the  Calcomp  computer  ploter  when  every  data  point  is  plotted. 

In  this  study,  the  planetary  bias  is  approximated  as  a  constant 
parameter  and  its  value  is  estimated.  But,  the  estimation  errors  and  the 
conditional  variances  remain  virtuall.y  constant  with  the  onboard  angle  mea¬ 
surements  during  the  time  period  of  interest.  Furthermore,  the  difference 
between  results  obtained  with  the  nonlinear  filters  and  the  EK  Filter  are 
negligible.  Consequently,  the  estimation  errors  and  conditional  variances 
corresponding  to  the  planetary  bias  are  not  presented  .in  this  report . 

Since  the  nonlinear  filters  are  compared  with  the  EK  filter,  the 
estimation  errors  and  conditional  variances  for  the  EK  filter  in  the  Simula¬ 
tions  1  through  8  should  be  identical.  However,  the  actual  figures  are  not 
identical  because  of  scale  factors. 


The  EK  filter  performs  adequately  up  to  three  days  and  thereafter 


*.iy 

becomes  unstable.  Actually,  right  after  the  three  day  period,  the  condition¬ 
al  variances  decrease  drastically  and  the  estimation  errors  take  several  sharp 
oscillatory  spikes  during  a  short  period  of  time  while  drifting  away  from 

zero.  The  sharp  decrease  in  conditional  variance  is  attributed  to  the  fact 
*  T 

that  h  Ph  dominates  R  in  Eq.  (1.35)  and,  hence,  the  negative  term  in 

X  X 

Eq.  (1.33)  will  be  quite  large.  The  position  and  velocity  estimation  errors 
and  the  conditional  variances  for  the  EK  Filter  are  shown  in  Figs.  8-a,  8-b, 
10-c,  and  10-d. 

Immediately  after  encounter  (12.7  days),  the  velocity  estimation 
errors  remain  at  a  relatively  constant  level  and,  hencr. ,  the  position  esti¬ 
mation  errors  grow  linearly,  in  an  unbounded  manner,  and  divergences  occur 
eventually.  The  velocity  estimation  errors  become  extremely  unstable  shortly 
after  encounter  and  the  magnitude  oscillates  several  times  with  ..iarp  spikes. 
This  phenomena  is  not  seen  in  the  figures  shown  here  because  of  the  fact 
that  every  third  data  point,  instead  of  every  data  point,  is  plotted.  In 
particular,  the  velocity  components  of  the  conditional  variances  are  very 
small  after  encounter,  and  the  filter  becomes  saturated.  Therefore,  the  ob¬ 
servations  taken  after  encounter  :annot  improve  the  estimate  very  much.  The 
characteristic  of  poor  estimation  after  encounter  is  an  indication  of  the 
importance  of  the  pre -encounter  navigation. 

Figs.  8-a  and  8-b  show  the  position  and  velocity  estimation  errors 
for  the  Simulation  1  which  compares  the  Filter  1  (or  MGSO  Filter)  and  the  EK 
Filter.  Both  Filters  perform  adequately  up  to  three  days  and  there  are  no 
significant  differences  between  them.  After  three  days,  the  EK  Filter  be¬ 
comes  unstable.  However,  the  MGSO  Filter  performs  properly  up  to  encounter. 
Both  Filters  diverge  after  encounter. 
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It  is  interesting  to  notice  that  the  M6S0  Filter  obtains  a  more 
accurate  estimate  of  the  Y  and  Z  components  and  less  accurate  estimate 
of  the  X  components  after  encounter  than  the  KK  Filter  does. 

Although  the  corresponding  conditional  variances  are  not  specifi¬ 
cally  shown  here,  they  are  almost  identical  to  the  ones  given  in  Figs.  11-d 
and  11-e.  It  is  interesting  to  note  that  the  conditional  variances  for  the 
MGSO  Filter  are  considerably  larger  than  those  of  the  EK  Filter  especially 
in  the  region  from  three  tc  thirteen  days,  during  which  the  MGSO  Filter  es¬ 
timates  surprisingly  better  than  the  EK  Filter. 

Figs.  9-a  and  S-b  show  the  position  and  velocity  estimation  errors 
for  Simulation  2  which  compares  Filter  2  and  the  EK  Filter.  Filter  2  per¬ 
forms  considerably  better  than  the  EK  Filter  throughout  the  entire  region. 

For  future  discussion  the  oscillations  around  encounter  are  emphasized  here. 
The  conditional  variances  are  identical  to  the  one  shown  in  Figs.  11-d  and 
11-e,  which  also  correspond  to  Simulations  1  and  u. 

From  Table  1,  it  can  be  seen  that  the  only  difference  between  Fil¬ 
ter  1  (or  MGSO  Filter)  and  Filter  2  is  that  the  dynamic  second  order  (DSO) 
term  is  dropped  in  Filter  2.  Therefore,  the  comparison  of  Filter  1  and  Fil¬ 
ter  2  shows  the  effect  of  DSO  term.  As  seen  in  Figs.  9-a,  8-!>,  9-a,  and  9-b, 
the  effect  of  DSO  term  has  c  significant  effect  after  encounter.  By  dropping 
the  DSO  term  from  the  MGSO  Filter  (Filter  t),  far  better  estimate  is  ob¬ 
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tained. 


A  number  of  simulations  indi.uie  that  the  DSO  term  is  very  sensi¬ 


tive  to  the  initial  covariance  matrix.  For  the  larger  values  of  the  initial 
variances,  less  satisfactory  estimates  are  obtained.  The  simulations  indi¬ 
cate  that  the  MGSO  Filter  diverges  while  th-_  LK  Filter  yields  «  convergent 
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estimate  whenever  large  initial  variances  are  used  with  relatively  large 
state  noise  Oq  ,  For  small  initial  variances  and  small  state  noise  Oq  , 
the  differences  between  the  two  Filters  MGSO  and  EK  Filters  are  negligible. 
This  implies  that  none  of  the  DSO,  0S0  and  KGC  terms  are  important.  Appa¬ 
rently,  most  of  the  orbit  determination  problems  which  are  not  influenced 
by  a  state  noise  u  fall  in  this  category  end  the  differences  between  the 
EK  Filter  and  the  MGSO  Filter  are  negligible.  However,  there  appears  to  be 
a  range  in  which  the  initial  variances  can  be  so  chosen  that  the  effect  of 
DSO  term  improves  the  filter  performance.  But,  it  may  not  be  easy  to  select 
such  an  initial  covariance  matrix  in  a  complex  multi-dimensional  problem, 
because  the  chosen  set  of  initial  covariance  matrix  may  very  well  cause  the 
DSO  term  to  affect  the  filter  in  such  a  way  that  the  estimate  of  some  com¬ 
ponents  can  be  improved  while  the  estimates  of  other  components  is  degraded. 
An  example  of  this  situation  is  shown  in  Figs,  8-a  and  8-b. 

Usually,  if  the  EK  Filter  converges,  i.e. ,  if  the  conditional  var¬ 
iances  remain  small,  the  MGSO  Filter  acts  like  the  EK  Filter  This  is  due 
to  the  fact  that  the  effect  of  the  DSO  term  cun  be  overridden  by  the  small 

j 

;  variances  associated  with  the  observations.  In  contrast,  if  the  covariance 

reduction  caused  by  the  observations  cannot  override  the  effect  of  the  DSO 
;  term,  which  will  occur  when  the  initial  variances  and  dynamic  state  noise 

are  large,  then  the  MGSO  Filter  diverges  because  of  the  DSO  term  even  though 
the  EK  Filter  converges.  A  large  conditional  variance  allows  the  estimate,  to 
depart  from  the  true  trajectory  because  o’"  the  DSO  term  and  cause  a  bad  pre¬ 
dicted  observation  and,  consequently,  large  observation  residual  which  will 
lead  to  filter  instability  and  divergence. 

Figs,  li-b  and  il-c  show  the  estimation  errors  for  Simulation 
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4  which  compares  Filter  4  with  the  EK  Filter.  As  seen  in  Table  1,  Filter  4 

includes  only  the  KGC  term  and  excludes  the  DSO  and  0S0  terms.  Filter  4  is 

referred  to  as  the  Kalman  Gain  Compensated  (KGC)  Filter. 

The  comparison  of  the  filter  performance  with  Filter  2  will  show 

clearly  the  effect  of  the  0S0  term  and  the  comparison  with  the  EK  Filter 

reflects  the  effects  of  the  KGC  term.  From  Figs.  11-b  and  ll-c,  it  is  seen 

that  the  Filter  4  (KGC  Filter)  estimates  show  excellent  agreement  with  the 

true  trajectory  throughout  the  entire  region.  Both  the  KGC  Filter  and  the 

» 

EK  Filter  appear  to  be  identical  for  the  first  three  days.  After  three  days, 
the  L'K  Filter  diverges.  Although  the  EK  Filter  performs  poorly  over  almost 
the  entire  region  except  for-  the  first  three  days ,  the  poor  performance  after- 
encounter  results  from  the  behavior  which  occurs  from  three  days  to  encounter. 
The  accumulated  large  estimation  errors  at  encounter  influence  the  estimate 
throughout  the  remainder  of  the  period. 

The  phenomena  above  can  be  explained  as  follov/s:  the  conditional 
variances  become  quite  small  after  encounter,  and  the  filters  become  insen¬ 
sitive  to  observations.  Therefore,  a  filter  that  can  estimate  accurately 
around  encounter  can  retain  an  accurate  estimate  after  encounter.  Similarly, 

any  filter  which  performs  in  an  unsatisfactory  manner  around  encounter  will 

« 

yield  an  inaccurate  estimate  after  encounter. 

Fig.  11-a  shows  the  observation  residuals  for  Simulation  4.  The 
observation  residual  pattern  for  the  EK  Filter  starts  to  grow  from  three 
days  and  is  influenced  by  a  large  spifce  around  encounter.  After  encounter 
the  residual  patterns  for  the  EK  Filter  and  the  KGC  Filter  remains  identical 
to  each  other. 


The  conditional  variances  are  shown  in  Figs.  13-d  and  tl-e.  It  is 
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interesting  to  note  that  the  conditional  variances  for  the  EK  Tilter  undergo 
a  sharp  reduction  at  a  time  of  about  four  days  following  which  the  estima¬ 
tion  errors  begin  to  drift  away.  The  conditional  variances  for  the  ’'GC  Filter 
retain  a  larger  value  for  the  period  of  time  from  four  to  twelve  days  during 
which  the  poor  performance  of  the  EK  Filter  has  been  accumulated. 

Fr-om  the  simulations ,  it  was  noticed  that  the  KGC  term  is  negli- 
gible  at  the  beginning  in  comparison  with  the  other  terms  h  £h  and  R  in 

X  X 

Eq.  (2.57).  But,  it  grows  rapidly  and  becomes  the  dominating  term  from  three 
to  twelve  days.  As  a  matter  of  fact,  the  maximum  value  of  the  KGC  term  is 
about  ten  times  larger  than  the  other  t \io  terms.  The  KGC  term  becomes  agai’* 
pjgligible  after  the  encounter.  The  above  fact  implies  that  the  observation 
nonlinearity  is  very  severe  from  three  to  twelve  days.  The  observations 
outside  this  region  appear  to  be  relatively  linear.  The  severe  observation 
nonlinearity  near  encounter  is  seemed  to  be  caused  by  tlv  dynamic  nonlinearity. 
The  same  investigation  was  made  on  the  other  type  oi  observations.  The  sun- 
planet  angle  measurement  has  almost  the  same  characteristics  as  the  star- 
planet  angle  measurement.  However,  the  range  and  range-rate  observations  do 
not.  appear  to  be  influenced  by  the  second  order  terms  and  the  KGC  terms  for 
both  observations  are  negligible.  Hence,  no  difference  between  the  KGC  Fil¬ 
ter  and  the  EK  Filter  is  seen. 

From  Eq.  (2.57),  it  can  be  seen  that  a  large  KGC  term  yields  a 
smaller  optimal  gain  K  than  that  which  results  in  the  EK  Filter.  Hfnce, 
a  smaller  conditional  covariance  reduction  occurs  and  a  larger  posterior  covar¬ 
iance  matrix  results,  as  can  be  seen  from  Eq.  (2.59).  In  Figs.  11-d  and  11-e 


exact!,  the  same  phenomena  described  above,  happens  in  the  region  from  tour 
days  to  encounter . 
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Between  observations ,  the  conditional  variances  vary  according  to 
'.ne  differential  Eq.  (5.55)  and  the  direction  of  change  depends  on  the  signs 
of  the  Jacobian  matrix  f  .  At  the  time  of  the  observation,  Eq.  (2.59) 
governs  the  conditional  covariance  matrix  reduction.  The  conditional  var¬ 
iances  increase  only  through  the  dynamics,  namely,  the  signs  of  f  and  de- 
by  either  the  dynamics  or  observations,  namely  h  .  For  example, 
the  reduction  for  the  EX  Filter  around  four  days  is  attributed  to  the  ob¬ 
servations  and  the  one  near  twelve  days  appears  to  be  due  to  the  sign  changes 

of  f  . 
x 

It  appears  that  near  encounter,  the  dynamic  nonlinearity  overrides 
the  information  gained  through  the  observations.  Physically,  this  means  that 
a  severe  dynamic  nonlinearity  causes  bad  predicted  observations  and  observa¬ 
tion  nonlinearity.  Therefore,  large  observation  residuals  are  inevitable. 

In  this  situation  one  can  follow  one  of  two  procedures: 

1.  Discard  the  observations  during  this  period. 

2.  Try  to  update  the  apriori  estimates  with  larger  gains,  K  . 

It  appears  that  the  EK  Filter  follows  the  second  course  while  the  KGC  Filter 
takes  the  first  course.  The  KGC  Filter  yields  a  large  value  for  the  condi¬ 
tional  covariance  matrix,  and  hence,  leads  to  a  small  value  for  the  gain  K 
because  of  the  KGC  term.  This  means  that  the  KGC  Filter  places  less  weight 
on  each  of  the  observations  obtained  during  the  period  of  time  in  which 
dynamic  nonlinearity  is  very  severe. 

In  addition  to  the  large  value  for  K  ,  the  observation  residual 
is  so  large  during  the  brief  period  of  time,  as  seen  in  Fig.  li-a,  that  the 
correction  to  the  apriori  estimate  in  the  EK  filter  becomes  e,  :essivsly 
large  and  a  poor  posterior  estimate  results.  The  conditional  variances 


for  llio  KGC  Filter  depend  more  on  the  state  dynamics  in  the  region  of  a 
severe  dynamic  nonlinearity.  The  sharp  decrease  of  the  conditional  variances 
near  encounter  is  not  due  to  the  observations  but  due  to  the  sign  changes 
of  the  Jacobian  matrix  f  .  h  large  integration  step  size  often  causes  a 
negative  variance  near  encounter.  The  sign  change  of  the  Jacobian  matrix 
f  incorporated  with  the  large  conditional  variances  yields  a  negative 

X 

slope  for  the  conditional  variances  which  can  result  in  negative  variances. 

The  conditional  variances  for  the  EK  Filter  depend  largely  on  the 
observation,  namely  hx  in  the  region  where  dynamic  nonlinearity  is  high. 

The  sharp  decrease  of  the  conditional  variances  for  the  EK  Filter  around  four 
days  is  an  indication  that  the  reduction  by  the  observation  overrides  the 
increase  due  to  the  state  dynamics.  Unusually  large  reductions  of  the  con¬ 
ditional  variance  which  occur  in  EK  Filter  in  the  early  stage  of  application, 
is  attributed  to  this  phenomena.  However,  it  can  be  a  nuisance  if  the  pos¬ 
terior  estimate  is  still  far  away  from  the  true  state  even  after  the  condi¬ 
tional  variances  reduced  to  a  zero  level.  The  most  interesting  observation 
is  that  the  EK  Filter  becomes  unstable  and  the  estimate  starts  to  drift  away 
from  the  true  trajectory  whenever  the  conditior.al  variances  are  reduced  sharp¬ 
ly.  Another  interesting  observation  is  that  *-ach  of  the  Filters  (1,  2  and  4) 
which  include  the  KGC  term  have  almost  identical  conditional  variances  snown 
in  Figs.  11-d  and  11-e,  and  have  very  similar  observation  residual  patterns, 
as  shown  in  Figs.  il-a. 

Figs.  12~a  and  12-b  show  that  the  estimation  errors  for  the  Simu¬ 
lation  5  which  reflects  the  characteristics  of  Filter  S.  Filter  5  5 ncl.udas 
only  the  QSO  term.  Filter  $  performs  better  than  the  EK  Filter  from  three  to 
nine  days.  Apparently,  the  observation  ..econd  order  improves  the  performance 


66 


of  Filter  5  by  using  the  proper  sign.  The  poor  performance  of  Filter  5  be¬ 
tween  nine  and  twelve  days  reflects  the  fact  that  the  0S0  term  influences 
the  filter  with  the  wrong  sign.  The  poor  estimates  of  Filter  5  after  en¬ 
counter  are  due  to  the  propagated  effect  of  the  poor  estimate  at  encounter. 

Figs.  10-a  and  10-b  show  the  estimation  errors  for  Simulation  3 
which  compares  Filter  3  with  the  EK  Filter.  Filter  3  contains  both  the  DSO 
and  the  0S0  terms;  By  comparing  Figs.  10-a  and  10-b  with  Figs.  12-a  and 
i2-b,  the  effect  of  the  DSO  term  is  shown  significantly  after  encounter. 

The  estimation  errors  for  Filter  3  are  considerably  larger  than  those  of 
Filter  5  after  encounter.  The  difference  would  be  the  negative  contribution 
of  the  DSO  term  in  the  Filter  3,  i.e.,  the  DSO  term  affects  the  filter  with 
the  wrong  sign.  The  conditional  variances  are  shewn  in  Figs.  10-c  and  10-d. 
These  figures  also  represent  the  conditional  variances  resulting  in  Simula¬ 
tion  5.  It  is  interesting  to  note  that  both  Filters  (3  and  5)  do  not  con¬ 
tain  the  KGC  term  and  the  conditional  variances  are  very  similar  to  the  one 
given  by  the  EK  Filter  and  quite  different  from  those  of  Filters  1,  2  and  4 
which  include  the  KGC  term.  Fig,  10-e  represents  the  observation  residuals 
of  Simulation  3.  This  observation  residual  pattern  which  contains  a  large 
spike  (even  larger  than  that  of  the  EK  Filter)  is  seen  also  in  Simulation  5. 
The  smooth  residual  pattern  (Fig.  li-a)  of  the  Filters  1,  2  and  '<  is  primar¬ 
ily  attributed  to  the  presence  of  the  KGC  term  anJ  its  side  effects. 

Simulations  0 ,  7  and  8  are  conducted  mainly  to  describe  the  char¬ 
acteristics  of  the  modified  truncated  s-cood  order  (MTSO)  filter  which  is 
designated  as  Filter  8.  rov  Simulations  6,  7  and  8,  only  the  X  components 
of  the  position  and  velocity  j:  *i mat ion  errors,  observation  res Final  and  the 

conditional  variances  V  and  V  are  shown  in  the  corresponding  figures. 

1 1 
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The  other  components  exhibit  similar  characteristics  and  are  no':  shown  in 
this  report  to  eliminate  unnecessary  space.  The  only  difference  between  the 
MGSO  Filter  and  the  MTSO  Filter  is  that  the  KGC  term  enters  with  a  plus  sign 
of  one-half  in  the  first  Filter  and  with  a  minus  sign  of  one-fourth  in  the  e 
second.  As  previously  pointed  out,  the  KGC  term  is  negligible  at  the  begin¬ 
ning  but  grows  rapidly  up  to  ten  times  the  value  of  the  combination  of  the 
'  T 

other  terms,  h  i!h  and  R  ,  as  shown  in  expression  (2.57).  Following  en- 
counter,  the  value  of  the  KGC  term  reduces  in  value.  From  the  above  charac¬ 
teristics  of  KGC  terms,  it  is  easy  to  conclude  that  the  MTGO  Filter  contains 
a  potential  signularitv.  The  optimal  gain  K  given  by  Eq.  (2.57)  with  a 
minus  sign  of  one-fourth  of  KGC  term  instead  of  plus  sign  approaches  plus' 

and  minus  infinity  as  the  KGC  term  approaches  the  sum  of  tne  ether  two  terms, 
*  T 

h  Ph  and  R  from  below  and  above.  In  addition,  when  the  cptiral  gain  is 
very  large  with  a  positive  sign,  the  posterior  conditional  variance  becomes 
negative.  The  phenomena  is  cl»arlv  'fleeted  in  Filters  6,  7  and  8,  and  is 
shown  in  Figs.  13-a,  13 -b,  14-a,  it.--  ,  15-a,  and  15-b. 

In  general,  tne  MGSO  I"  ' ter  keeps  the  ''onditional  variances  larger 
than  those  of  the  EK  Filter.  The  MT;  •  Filter,  in  contrast,  has  a  tendency 
to  keep  the  conditional  va  -'ances  u ,  IlL  ^  i,  an  rho' *  of  tie*  EK  Filter.  How¬ 
ever,  the  variance  b- cor.--'''  meaning! .  as  tne  term  drive,  the  variance 
to  a  negative  quantity  in  Hirers  . ,  «nd 
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Figure  9-a.  Position  Estimation  Errors  for  the  Simulation 
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Figure  14-b.  Conditional  Variances  (V.^V..)  for  the  Simulation 
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4.4  Applications  of  the  KGC  Filter  to  the  Hypobolic  Orbit 

In  the  previous  section,  several  filters  were  discussed  and  che  DSO 
and  0S0  terms  were  characterized  in  conjunction  with  the  KGC  term.  Since 
the  KGC  term  is  the  only  term  which  improves  the  fi-ltc:  performance*  the  KGC 
filter  which  includes  only  the  KGC  term  is  designated  as  the  best  filter 
among  those  listed  in  Table  1.  > 

The  KGC  filter  is  further  tested  through  numerous  simulations.  The 
Simulations. 9;  through  14  are  conducted  on  a  hypobolic  orbit  with  the  sun- 
planet  angle  measurements .  The  simulations  are  designated  to  determine  the 
effects,  of  state  noise  covariance  matrix  Q  ,  initial  state  errors  > 

-integration  step  size.,  observation  rate.,  initial  covariance  V0/0  and  ob¬ 
servation  noise  covariance  matrix  R  .  Simulation  10  is  the  reference  case 

to  which  all  other  simulations  are  compared.  The  input  data  are  given  in 

-  . 

Tables  ’?■  and  8. 

Simulation  9  is  specifically  designed  to  illustrate  the  effect  of 
'using  two  different  Q's  in  Eq.  (2.55).  o^T  is  the  true  standard  deviation 
of  the  state  noise  u  which  is  used  to  compute  the  true  trajectory,  is 

the  apriori  standard  deviation.  The  square  of  is  used  in  Eq.  (2.55) 

£or-  the  estimation  procedure.  It  is  common  practice  to  inc3.ude  a  Q  in  Eq-. 
(2. 55) ,  although  there  may  be  no  state  noise  u  assumed.  This  procedure  is 
followed  to  keep  the  value  of  the  conditional  variances  above  a  certain  level 
so  that  the  filter  can  maintain  a  reasonable  gain  K  and,  hence,  will  be 
sensitive  to  the  observations.  Usually  the  EK  filter  reduces  the  conditional 
covariance  very  rapidly  after  a  few  observations  ere  made,  and,  hence,  the 
filter  becomes  saturated  and  insensitive  to  the  observatiens. 

Figs.  16-a  and  16-b  show  the  position  and  velocity  estimation 
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errors,  respectively,  for  Simulation  9.  The  conditional  variances  and  ob¬ 
servation  residual  are  almost  identical  to  the  ones  gi^en  in  Figs.  17-c, 

17-d  and  17-e  which  are  obtained  by  using  the  same  and  o  .  The 

observation  residual  pattern  should  not  be  assumed  to  be  zero  except  during 
the  brief  time  interval  in  which  the  spike  occurs.  The  non-zero  value  of 
the  residuals  do  not  show  up  on  the  scale  used  to  plot  the  results. 

Since  the  sun-plahet  angle  measurement  is  restricted  to  the  eclip¬ 
tic  plane,  the  measurement  does  not  include  very  much  information:  on  the  Z 
components  of  position  and  velocity.  This  fact  is  reflected  in  the  figures 
related  to  the  Z  components  from  zero  to  ten  days.,-  Both  the  EK  filter 
and  the  KGC  filter  perform  reasonably  well  up  to  eight  days.  As  a  matter 
of  fact,  they  are  almost  identical.  The  EK  filter  starts  to  drift  away  from 
the  true  trajectory  after  eight  days.  Around  twelve  days,  the  EK  filter  be¬ 
comes  extremely  unstable,  oscillates  several  times  with  sharp  spikes  during 
this  short  period  of  time  and  then  diverges  eventually.  The  oscillatory 
spikes  near  encounter  are  not  shown  in  the  figures  simply  because  only  every 
third  data  point  is  shown. in  the  figures.  As  seen  in  Fig.  3.7-e,  the  esti¬ 
mate  with  the  EK  filter  is  influenced  by  a  spike  in  the  observation  residual 
pattern  around  encounter.  The  actual  observations  which  depend  mainly  oh 
the  true  states  and  the  small  observation  noise  contain  equally  good  infor¬ 
mations  at  any  time.  But  the  predicted  observation  depending  on  the  current 
estimate  can  be  quite  erroneous.  The  erroneous  observation  residuals  around 
encounter  are  incorporated  with  a  large  optimal  gain,  K  ,  during  the  same* 
period  and  the  update  of  the  apriori  estimates  is,  consequently,  too  large, 
causing  the  EK  filter  to  diverge. 

The  KGC  filter  keeps  the  conditional  variances  large  in  the  critical 
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period  of  time  from  eight  to  twelve  days.  As  a  result,  the  KGC  term  dominates 

a  ij> 

the  other  two  terms,  h  Ph  and  R  ,  and  the  values  of  the  optimal  gain,  K  , 
remain  small.  Through  this  effect,  the  KGC  term  makes  the  filter  insensitive 
to  the  observations  which  could  be  erroneous ,  in  the  region  where  the  dynamic 
nonlinearity  is  very  severe. 

Simulation  10  shows  the  effect' of  using  the'same  values'  for 
and  Oqa  .  it  is  noticeable  that  Simulation  9  wnich  user  two  different  val¬ 
ues  for  anr  and  on,  ,  yields  slightly  smaj  estimation  errors  for  both 
x*  QA 

the  EK  and  KGC  filters  than  those  of  the  Simulation  10  which  uses  the  same 
value.  Since  the  sun-planet  angle  measurement  is  restricted  to  -the  ecliptic 
plane,  the  EK  filter  experiences  severe  nonlinearity  effects  on  the  Z  com¬ 
ponents.  The  effect  can  be  seen  in  the  Z  components  of  the  conditional 
variances  shown  in  Figs.,  ,17 -c  and  17-d. 

Simulation  11  shows  the  effect  of  the  initial  state  errors.  The 
errors  are  chosen  ten  times  larger  than  those  of  Simulation  10.  The  conver¬ 
gence  characteristics ,  except  the  Z  components ,  appear  to  be  well  behaved 
for  both  filters  immediately  after  taking  the  observations.  However,  the 
EK  filter  displays  instability  and  divergence  characteristics  around  en¬ 
counter,  although  the  KGC  filter  continues  with  an  accurate  estimate  through¬ 
out  this  extremely  nonlinear  region. 

Overall,  the  estimation  error  patterns  shown  in  Figs.  18-b  and  18-c 
are  very  much  the  same  as  those  of  Simulations  9  and  10,  except  during  the 
first  few  days.  The  observation  residual  pattern  is  shown  in  Fig.  18-a,  and 
the  conditional  variances  are  almost  identical  to  the  one  given  in  Figs,.  17-c 
and  17-d.  Again,  the  residuals  outside  the  spike  zone  do  not  show  up  because 
of  the  relatively  small  size  compared  with  the  spike. 
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Since  the  nominal  trajectory  changes  very  rapidly  around  encounter 
in  the  hypo'oolic  orbit  and  near  periapsis  in  the  elliptic  orbit,  the  dynamic 
nonlinearity  appears  to  be  very  severe  in  these  regions  and,  consequently, 
so  is  the  observation  nonlinearity.  In  order  to  minimize  the  effects  of  the 
nonlinearities,  a  variable  integration  step  size  and  observation  rate  are 
adopted  in  the  Simulation  12.  Initially,  the  integration  step  size  and  the 
observation  intervals  are  set  to  1/10  day,  and  then  the  actual  integration 

step  size,  AT.  ,  and  observation  interval,  AT  ,  are  determined  as  follows: 

.  '*■  -  '-ow-  - 

AT.  =  AT  =  (1710)  *  Integer  Value  of  (r/r  )  ,  r  >  r 

i  o  o  -  o 

or  -  -  -  .  j 

AT...  =  AT.  =  (1/10)  /  Integer  Value  of  (r  /r)  ,  f  >  r 

"1  .0  o  o 

where  r^  is  the  initial  distance  between  the  spacecraft  and  the  target 
planet,  Jupiter,  and  r  is  the  current  distance.  As  seen  in  Fig.  7,  AT,, 
arid  ATq  ,  for  example,  become  nearly  1/20  of  the  initial  step  size,  1/10 
day  hear  encounter  in  the  hypobolic  orbit. 

Figs.  19-a  and  19-b  shov?  the  estimation  errors  for  the  position 
arid  velocity  components.  An  interesting  fact  about  the  error  pattern  of 
Simulation  12  is  -chat  the  signs  are  reversed  when  compared  with  those  of 
Simulation  10.  Figs.  19-c  and  19-d  show  the  conditional  variances  and  a 
slightly  different  pattern  is  seen  near  encounter  when  compared  with  those 
in  Figs*  17-c  and  17-d. 

Simulation  13  is  conducted  with  a  ten-times  larger  initial  condi¬ 
tional  variances  of  the  velocity  components  than  those  of  Simulation  io.  Due 
to  the  larger  initial  conditional  variances ,  the  KGC  term  starts  influencing 
the  KGC  filter  earlier  than  it  does  in  Simulation  10,  Figs.  20-c,  20-d  and 


20-e  represent  the  conditional  variances  and  observation  residuals,  respect¬ 
ively.  Estimation  errors  are  shown  in  Figs.  20-a  and  20-b. 

Simulation  14  shows  the  effect  of  a  large  observation  noise  stand¬ 
ard  deviation  oR  in  the  KGC  filter.  A  value  ten  times  larger  than  the 
value  of  oR.  which  was  used  in  Simulation  10  was  adopted  as  an  observation 
noise  standard  deviation  in  Simulation  14. 

r 

It  is  found  in  the  EK  filter  which  docs  not  include  the  KGC  term 

.  T  / 

that  hxPh^  dominates  R  in  Eq.  ( l|.  35)  during  the  early  stage  of  estitna- 
tion.  Later  the  values  of  R  dominate.  From  the  above  observation,  it  is 

understood  that  the  KGC  term  which  is.  negligible  when  compared  with  the 

„  T  _ 

other  two  terms,  h^Ph,  and  R  in  Eq.  (2*57),  cannot  affect  the  performance 

of  the  KGC  filter  very  much.  But  in  the  region  where  the  KGC  term  dominates 

a  t 

the  other  two  terms  and  R  is  larger  than  hvPh  ,  the  effect  of  a  large 
R  bhows  up.  For  example,  consider  the  period  from  nine  to  thirteen  days 
in. Figs.  21-b  and  21-c.  After  encounter,  the  value  of  the  KGC  term  diminishes 
due  to  the  combined  effect  of  small  conditional  variances  and  observation 
second  partials  h  and  thereafter  the  KGC  term  has  virtually  no  influence 

XX 

on  the  KGC  filter.  Therefore,  the  KGC  filter  performs  like  the  EK  filter 
after  thirteen  days.  For  a  small  oR  ,  the  effect  of  the  KGC  term  becomes 
very  significant  and  the  KGC  filter  is  very  desirable  whenever  a.  smaller  oR  , 
or  equivalently  accurate  observation,  is  available. 

The  corresponding  conditional,  variances  and  observation  residuals 
are  shown  in  Figs.  21-d,  21-e  and  21-a,  respectively. 
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Figure  19-b.  Velocity  Estimation  Errors  for  the  Simulation  12 
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*1;  !>  Appi.1  frit  Ion  of  the  KGC  Filter- -to  the  Ullipl  ic  Orbit  _  , 

The  KGC  Filter  is  applied  to  the  problem  of  determining  an  ellip¬ 
tic  orbit  around  Jupiter  with  sun-planet,  star-planet  and  range-rate  plus 
sun-planet  ..angle  measurements .  Since  the  variation  of  the  Z  components  in 
the  nominal  trajectory  shovm  in  Fig.  6  is  very  small  compared  with  the  X 
and  Y  components,  the  initial  state  errors  for  the  Z  components  of  posi¬ 
tion  and  velocity  are  chosen,  as  one-tenth  of_-those  of  X  and"Y“  "components,, 
i.e. ,  X  and  Y  components  of  position  and  velocity  estimation  errors  are  , 
initially  chosen  as  103  km  and  1CT4  km/sec  ,  respectively,  and'  Z  com¬ 
ponents  as  102  km  and  10" 5  km/sec  .respectively. 

Figs.  22-a  and  22-b  show  the  estimation  errors  for  Simulation  15. 
Unlike  the  hypobolic  orbit,  the  EK  Filter  gradually  drifts  away  after  fifteen 
days  during  which  time  the  dynamic  nonlinearity  affects  grow  large.  The  con¬ 
ditional  variances  are  shown  in  Figs.  22-c  and  22-d.  Fig.  22-e  shows  the 
observation  residual  pattern  which  consists  of  every  tenth  data  point  con¬ 
nected  with  straight  lines.  Simulation  18  is  designed  to  see  the  performance 
of  the  KGC  Filter  for  a  long  period  of  time.  The  period  of  estimation  is 
extended  to  62  days  which  is  more  than  two  revolutions  of  the  elliptic  orbit. 
The  same  input  data  as  that  used  in  Simulation  15  are  used  for  Simulation  18. 
The  estimation  errors  for  the  position  and  velocity  are  shown  in  Figs.  25-a 
and  25-b.  For  the  first  30  days,  the  estimation  error  patterns  of  Figs.  25-a 
and  25-b  match  identically  with  the  ones  shown  in  Figs.  22-a  and  22-b.  An 
interesting  fact  is  that  both  the  EK  Filter  and  the  KGC  Filter  exhibit  a  per¬ 
iodicity  in  the  estimation  errors.  However,  the  errors  for  the  EK  Filter 
grow  larger  during  the  second  revolution  and  reach  an  unacceptable  value. 
However,  the  KGC  Filter  realizes  a  very  accurate  estimate  throughout  the  en¬ 
tire  period  and,  as  a  matter  of.  fact,  the  estimation  errors  are  considerably 


smaller  during  the  second  revolution  than  those  during  the  first  revolution. 
Around  the  second  periapsis,  the  EK  Filter  estimation  errors  exhibit  spikes 
in  the  velocity  estimation  error  components  and  the  numerical  value  of  the 
.position  changes  sign.  The  conditional  variances  are  shown  in  Figs;  25-c  and 
25-d. 

From  the  examination  of  Figs.  25-a,  25-b ,  2b~c  and  25-d,  it  can  be 
seen  that  the  KGC  Filter  is -superior  to  the  EK  Filter  in  the  region  where 
the  conditional  variances  of  .the  KGC  Filter  are  larger  than  those  of  the  EK 
.  -Filter.  The  improvement  achieved  in  the  KGC  Filter  is  strictly  due  to  the 
effect  of  the  KGC  term  in  the  optimal  gain  K.  Fig.  25-e  shows  the  observa- 
*  tion  residual.  It  shows  a  couple  of  spikes  around  the  second  periapsis. 

For  the  blown  up  scale ,  the  residual  pattern  of  every  tenth  data  point  for 
the  first  30  days  matches- with  the  one  -shown  in  Fig.  22-e.i 

Simulation  16  is  conducted  with  sun-planet  angle  measurements .  The 
estimation  errors  shown  in  Figs.  23-b  and  23-c  reflect  the  characteristics 
of  the  sun-planet  angle  measurement.  Since  the  sun,  Jupiter  and  the  space¬ 
craft  are  all  almost  on'  the  ecliptic  plane,  the.  information  about  the  Z 
,  components  is  poor  and  the  EK  Filter  determines  a  poor  estimate  of  the  Z 
components  of  position  and  velocity.  Figs.  23-d,  23-e  and  23-a  show  the 
conditional  variances  ana  the  observation  residual  pattern ,  respectively. 

.  The  residual  pattern  is  obtained,  by  connecting  every  tenth  data  point  with  a 
straight  line .  .  r 

",  .  The  Simulation  17  is  conducted  with  two  kinds  of  observations,  i.e. 

rangerrate  plus  sun-planet  angle  measurements.  Figs.  24-a  and  24-b  show  the 
estimation  errors  and  the  conditional  variances  are  given  in  Figs.  24-d  and 
24-e.  The  Z  components  are  estimated  poorly  again.  Fig.  24~e  -shows  the 


two  observation  residuals:,  one  for  the  range-rate  and  the  other  for  the  sun- 
planet  angle  measurement.  With  the  data  of  Simulations  16  and  17,  the  filters 
were  examined  for  62  days.  Th  same  characteristics  discussed  in  Simulation 
18  can  be  found.  The  conditional  variances  vary  periodically  and  large  spikes 
can  be  found  in  t  je  observation  residual  and  velocity  estimation  errors  for 
the  EK  Filter.  The  EK  Filter  estimation  errors  are  incomparable  at  the  sec¬ 
ond  revolution.  Actually,  they  diverge  after  the  second  periapsis.  However, 
the  KGC  Filter  performs  exceptionally  well  through  out  the  entire  period  and 
even  better  at  the  second  revolution. 

Simulations  1  through  8  are  re-examined-  by  using  the  variable  ir.te- 
gration  step  size  and  observation  rate  discussed  in  Section  4,4.  The  char¬ 
acteristics  discussed  in  the  Section  4.3  are:  unchanged. 

It  is  found  that  variable  integration  step  size  and  rapid  observation 
rate  do  not  change  the  characteristics  of  Simulations  15,  16,  17  and  18  ex¬ 
cept  that  over  $11  estimation  errors  are  smaller-  than  those  of  the  constant 
integration  step  size  and  observation  rate.  However,  the  improvement  is  not 
significant  considering  the  computer  time.  Usually,  for  this  case,  the  EK 
Filter  reduces  the  estimation  errors  more  significantly  than  the  KGC  Filter 
does,  but  the  EK  Filter  solution  will  diverge  still. 
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*Estimation  performed  for  62  days  which  is  about  two  revolutions. 
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Figure  22-b.  Velocity  Estimation  Errors  for  the  Simulation  15 
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Conditional  Variances  of  Position  Estimation  Errors  for  the  Simulation  17 
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Figure  24-e.  Observation  Residuals  for  the-f  Simulation  17 
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Figure  25-b.  Velocity  Estimation  Error  for  the^ Simulation  18 
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CONCLUSIONS  AND  RECOMMENDATIONS 

5.1  Summary 

In  the  investigation  described  in  the  previous  discussion,  a  compar¬ 
ative  study  has  been  made  of  nonlinear  estimation  algorithms  and  their  appli¬ 
cation  to  the  orbit  determination  problem  for  interplanetary  spacecraft.  By 
using  the  properties  of  a  Martingale  series  and  Loeve’s  smoothing  properties, 
a  second  order  nonlinear  estimation  algorithm  is  derived.  The  algorithm  is 
shown  to  be  of  the  Gaussian  second-order  class  as  distinguished  from  the 
truncated  second-order  class.  Both  classes  of  second  order  filters  retain  a 
second  order  term  in  the  state  dynamics,  the  observation  state  relation  and 
in  the  optimal  weighting  matrix  (Kalman  gain),  respectively.  The  merits  of 
each  of  the  algorithms  as  well  as  the  influence  of  each  second  order  term  is 
evaluated  by  a  numerical  simulation  of  the  orbit  determination  for  a  Jupiter 
fly-by  and  Jupiter  orb iter  missions. 

5.2  Conclusions 

Based  on  the  results  of  extensive  numerical  simulations  on  the 
Jupiter  fly-by  and  Jupiter  orbiter  missions,  the  following  conclusions  can  be 
drawn  for  the  class  of  problems  considered  here: 

1.  The  system  dynamic  influences  the  performance  of  the  EK  Filter 

through 

i)  initial  conditions 

ii)  conditional  covariance  matrices 

iii)  predicted  observations 
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The  effect  of  the  initial  conditions  is  considered  to  effect  the  system  dy¬ 
namics  directly,  i.n  contrast  to  the  effects  of  the  other  two  factors,  which 
are  regarded  as  an  indirect  effect  on  the  performance  of  the  filter.  1't  is 
concluded  that  the  indirect  effects  of  the  system  dynamics  are  more  severe 
than  the  direct  effects,  especially  when  the  system  dynamics  are  highly  non¬ 
linear  . 

2.  The  effect  of  the  dynamic  second  order  (DSO)  term  cannot  be 
isolated.  If  the  conditional  variances  are  large ,  which  means  that  the  ex¬ 
tended  Kalman  Filter  does  not  perform  adequately,  the  effect  of  DSO  term  is 
very  severe  and  causes  the  second  order  filters  (MGSO  and  MTSO)  to  c:*'jrge, 
in  a  situation  when  the  EK  Filter  performs  reasonably  well.  In  contrast,  the 
small  conditional  variances,  which  imply  that  the  EK  Filter  works  very  well, 
do  not  reveal  any  differences  between  the  second  order  filters  (MGSO  and  MTSO) 
and  the  EK.  Filter.  \s  a  matter  of  fact,  there  is  no  reason  for  using  a  sec¬ 
ond  order  filter  if  the  conditiona  .ariances  are  small  and  tne  EK  Filter 
performs  adequately.  By  including  the  dynamic  second  order  (DSO)  term,  ap¬ 
proximately  30%  more  computer  time  is  required  than  that  required  by  the  EK 
Filter.  Furthermore,  any  filter  including  the  DSO  term  is  very  sensitive  to 
the  initial  covariance  matrix,  if  a  dynamic  nonlinearity  is  significant  at 
the  beginning.  It  is  interesting  to  note  that  Athans  et  al.  (27),  based  on 

a  different  example  problem,  concluded  that  the  DSO  term  is  the  major  fac¬ 
tor  in  improving  the  performance  of  the  MGSO  Filter,  Hence,  the  conclusions 
reached  in  this  investigation  regarding  the  DSO  term  should  be  regarded  as 
problem  dependent. 

3.  The  observation  nonlinearity  depends  on  the  type  of  observa¬ 
tions  and  the  dynamic  nonlinearity.  The  range  and  range-rate  observation  in 
the  Jupiter  fly-by  problem  are  regarded  as  relatively  linear  observations. 
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Like  the  DSO  term,  the  effect  of  the  observation  second  order  (0S0)  term  can¬ 
not  be  determined  for  all  problems.  The  same  conclusions  as  those  reached 
for  the  DSO  term  can  be  stated  for  the  0S0  term.  But  the  effect  of  the  OSO 
term  is  not  as  severe  as  the  effect  of  the  DSO  term  for  the  MGSO  Filter. 

4.  The  Kalman  Gain  Compensated  (KGC)  Filter  appears  to  give  a 
quite  acceptable  behavior  based  on  the  following  observations: 

i)  In  the  region  where  dynamic  nonlinearity  is  not  significant, 
the  KGC  Fill aFT.au ts  liVTTihe  extended  Kalman  filter, 

ii)  In  the  region  where  dynamic  nonlinearity  is  very  severe  and, 
consequently,  the  observation  residuals  are  large,  the  KGC 
Filter  down  weights  the  large  observation  residuals, 

iii)  The  effect  of  the  KGC  Filter  becomes  more  significant  when  the 
observation  noise  R  is  small  and  the  state  noise  Q  is  large. 
This  fact  implies  that  the  KGC  Filter  is  more  desirable  when 
the  observations  are  measured  accurately  and  when  the  dynamic 
noise  is  large  or  equivalently  when  the  dynamic  process  is 
poorly  modeled. 

iv)  The  KGC  Filter  is  very  stable  and  insensitive  to  the  dynamic 
nonlinearity  as  compared  with  the  EK  Filter. 

v)  The  KGC  Filter  maintains  an  accurate  estimate  for  the  highly 
nonlinear  type  of  observations  while  acting  like  the  extended 
Kalman  Filter  for  the  relatively  linear  type  of  observations, 

vi)  In  contrast  to  any  other  second  otder  filter,  implementation 
of  the  KGC  Filter  is  as  feasible  for  complex  problems  as  the 
extended  Kalman  Filter  is. 

5.  If  no  state  noise  is  assumed,  the  EK  Filter  works  adequately 
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and  the  differences  between  the  two  filters*  (EK  and  K6C)  are  negligible.  The 
absence  of  state  noise  implies  that  the  conditional  variances  cannot  be  too 
large . 

7.  for  the  Jupiter  fly-by,  the  extended  Kalman  Tilter  determines 
an  etdequate  estimate  up  to  a  period  of  encounter  minus  three  days.  However, 
the  estimate  diverges  around  encounter,  when  the  dynamic  state  noise  is  in¬ 
cluded  while  the  KGC  Filter  yields  accurate  estimates., 

8.  For  the  Jupiter  orbiter,  the  estimate  of  the  extended  Kalman 
Filter  drifts  away  gradually  from  the  true  trajectory  and  diverges  at  the 
second  revolution,  when  the  dynamic  state  noise  is  included  while  the  KGC 
Filter  yields  an  accurate  estimate. 

5.3  Recommendations  for  Future  Study 

The  research  reported  here  .is  an  indication  of  a  successful  appli¬ 
cation  of  an  approximate  nonlinear  filter  and  indicates  the  possibility  that 
the  Kalman  Gain  Compensated  Filter  can  be  applied  to  other  problems.  The 
following  studies  are  recommended: 

1.  Application  of  the  KGC  Filter  to  the  orbit  determination  prob¬ 
lems  associated  with  re-entry,  near-Earth  and  lunar  satellites.  Mariner  and 
Viking  missions,  should  be  carried  out.  In  particular,  application  of  the 
KGC  Filter  to  the  re-entry  and  ascent  phases  of  the  shuttle  navigation  prob¬ 
lem  is  recommended. 

2.  The  applicability  of  the  square  root  covariance  and* the  consider 
filter  versions  of  the  Kalman  Gain  Compensated  Filter  should  be  developed. 

3.  An  extended  study  of  nonlinear  estimation  algorithms  and  their 
applicability  to  the  orbit  determination  problems  should  be  made.  The  com¬ 
parative  study  should  be  made  in  the  frame-work  of  the  applicability  of  the 


methods  to  anticipated  orbit  determination  problems.  In  the  study,  particu¬ 
lar  attention .should  be  given  to  those  data  types  which  undergo  significant 
geometrical  changes  during  a  mission.  The  objective  of  such  a  study  would 
be  to  define  particular  missions  and  data  types  for  which  nonlinear  orbit 
determination  algorithms  will  yield  a  significant  improvement  over  the  es¬ 
timate  obtained  with  the  extended  Kalman  Filter. 

4.  Further  study  of  the  effects  of  the  dynamic  second  order  (DSO) 
term  and  the  observation  second  order  (0S0)  term  should  be  made. 

5.  An  .extensive  study  of  the  Gaussian  second  order  filter  and  the 
truncated  second  order  filter  should  be  made  in  the  direction  of  determining 
the  characteristics  of  the  random  forcing  term  in  the  covariance  equation. 


ArPENDIX  A 


THE  PARTIAL  DERIVATIVES 


The  first  order  partial  derivatives  f  and  h  and  the  second 
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order  partial  derivatives  f  and  h^x  which  appear  in  the  nonlinear  es¬ 


timation  algorithm  are  defined  in  this  appendix. 


From  Eq.  (3.5),  the  (9*9)  matrix  f  can  be  partitioned  as  follows 


4 

I 

* 

f  = 

X 

A21 

* 

A23 

4> 

- 1 

where  <i>  is  the  3^3  null  matrix  and  I  is  the  3x3  identity.  The  symmetric 


submatrices  and  A23  are  defined  as  follows: 
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The  second  order  partial  derivatives  are  defined  as  a  (3*9*9)  three- 


dimensional  array  and  each  layer  is  defined  as  a  symmetric  (9*3)  matrix.  The 


three-dimensional  array  can  be  pictured  as  shown  in  Fig.  26. 


Each  layer  is  further  partitioned  into  (3*3)  matrices  and  they  are  as  follows 
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where  <J>  is  a  (3*3)  matrix  with  identically  zero  elements.  Tha  sub-r.ati ices 


F^i  ,  Fj3  ,  and  F^  are  redefined  as  follows; 
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The  second  order  partial  derivative  h  is  a  (4x9*9)  three-dimensional 
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array  shown  in  Fig. '27  and  each  layer  is  a  (9x9)  symmetric  matrix. 

Each  layer  is  partitioned  into  (3*3)  matrices  and  they  are  as  follows: 


All  submatrices  with  the  exception  of  D  are  (3*3)  symmetric  matrices  and 
they  are  defined  as  follows: 
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APPENDIX  B 


CONDITIONAL  EXPECTATION  AND  SMOOTHING  PROPERTY 
Loeve  (34)  introduced  the  concept  of  "conditioning"  in  terms  of  sub 
o-fields  of  events,  and  further  conditional  probabilities  of  events  and  con¬ 
ditional  expectations  of  random  variables  "given  a  o-field  8"  as  8-measur- 
able  functions  defined  up  to  an  equivalence.  The  conditional  probability  of 
an  event  A  "given  an  event  B"  corresponds  to  that  of  the  frequencies  of  the 
occurrence  of  A  in  the  repeat  trials  where  B  occurs.  For'  every  event  A,  the 
relation 


PB  ’  P_A  =  PAB 

D 


(B-l) 


defines  the  conditional  probability  P(A/B)  of  A  given  B  as  the  ratio 
PAB/PB  ,  provided  B  is  a  nonnull  event.  In  a  more  mathematical  .-ox-ui,  the 
function  Pn  on  the  o-field  A  of  events ,  whose  values  are  P0/>  ,  AeA  , 
is  called  the  conditional  probability  of  A  given  B  .  3ince  P  on  A  is 
normed,  nonnegative  and  o-additive,  so  is  PD  on  A  ;  and  P  satisfies  the 
following  condition 

V2  =  1  »  PB  2  0  ’  Vi  =  >:PBAi 

Thus,  the  conditioning  expressed  by  "given  B"  means  that  the  initial  proba¬ 
bility  space  (n,A,P)  is  replaced  by  the  probability  space  (fl,A,P0)  .  The 
expectation  of  a  random  variable  X  on  this  new  probability  space  is  called 
the  conditional  expectation  given  B  and  is  denoted  by 


E[X/B]  = 


XdPr 
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mpb  *  L  mpb 


(B-2) 


whvro  B  is  the  complement  of  B  .  Knowing  that  =  0  on  {AB^,Ar.A} 


and  that  r 


1 


r—  !'  on  IAB,AcA)  from  Eq.  B-,t,  Eq.  B-?  becomes 


E(X/B]  = 


PB  J 


XdP 
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(B-3) 


B 


which  is  the  definition  of  the  conditional  expectation  of  X  given  B  .  The 
conditional  expectation  acquires  its  full  meaning  when  interpreted  as  values 


oi  functions  as  follows:  the  number  C[X/B]  is  no  longer  assigned  to  B  but 


to  every  point  of  B  ,  and  similarly  for  E[X/B  j  ,  so  that  we  have  a  two- 


C  C 

valued  function  on  f l  with  values  E[X/B]  for  weB  and  E[X/B  ]  for  weB 


More  generally,  let  { B_. }  be  a  countable  partition  of  ft  and  let  8  be  the 


minimum  o-field  over  this  partition.  Let  E  be  the  family  of  all  random  var¬ 
iables  X  whose  expectation,  E(.X]  ,  exists  so  that  their  indefinite  integrals, 
hence  conditional  expectations  given  any  nonnull  event,  exist.  Then  the  condi¬ 
tional  expectation  of  X  given  6  is  defined  as  the  following  elementary 
functions  (see  page  64  of  Loeve  for  the  definition)  up  to  an  equivalence 


E [ X/B 1 
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irf 

n  '  V 


i  XdP)l_  j  XeE 

PBj  V  Bj 


(B-4) 


The  above  is  the  constructive  definition  and  is  different  from  (B-3)  in  the 
sense  that  conditioning  is  given  as  a  o-field,  3,  instead  of  an  event  B  , 

Bc6  .  It  can  be  easily  seen  further  that  the  conditioning  car.  be  either  as  a 
random  variable  which  is  8-measurable  function  or  as  an  output  of  the  random 


variable.  If  the  partition  {B^}  is  not  countable,  the  above  constructive 


definition  is  not  applicable  and  rather  powerful  tool,  namely,  the  Radon- 
Nikodym  theorem  is  employed  and  the  descriptive  definition  is  followed,  let 
Pg  be  the  restriction  of  P  to  8  ,  defined  by 


P8B  =  PB  , 


Be6 


(B-5 ) 


then  the  conditional  expectation  E[X/6]  of  X  given  8  is  anv  S-measur- 


abel  function  whose  indefinite  integral  with  respect  to  P^  is  the  restriction 


m 


m 


~>,V,  , ^*3  j  JH  It  ■  . 


i  r>!> 

to  8  ol  the  indefinite  integral.  of  X  with  respect  to  1*  .  This  definition 
means  precisely  that,  for  every  XcX  ,  E[X/B]  is  defined  by 

E(X/8]dPp  =  XdP  ,  BeB  (B-6) 

B  -*8 

up  to  an  equivalence.  Loosely  speaking,  the  operation  E[X/B]  is  a  B-smooth- 
ing  and  some  of  its  important  properties  are  quoted  from  Loeve  without  proof. 

1.  On  every  nonnull  atom*  BeS  ,  C[X/8]  is  constant  and  its  value  E[X/B] 

is  the  average  of  the  values  of  X  on  B  with  respect  to  P  . 

2.  For  every  8  independent  of  the  o-field  8^  of  events  induced  by  X 

E[X/B]  =  EX  a.s. 

3.  Conditional  expectation  operator  E[*/8]  and  g-measurablc  factors  com¬ 
mute,  that  is,  if  X  is  8-:neasurable ,  then 

E[XY/8)  =  XE[Y/8]  a.s. 

4.  If  BcB'  ,  then 

E[E[X/S* ]/B]  =  E[X/8]  =  E[E[X/B]/B’ ]  a.s. 

It  is  interesting  to  note  that  for  the  "least  fine"  or  "smallest"  of  all  pos¬ 
sible  e-fields  BcA  ,  that  is,  for  8q  =  {^>,51}  ,  E tX/8Q ]  =  E[X]  almost 
surely,  which  means  that  unconditional  expectation  is  a  special  case  of  con¬ 
ditional  expectation  whose  conditioning  is  merely  the  least  fine  o-field. 

It  is  noted  also  that  any  deterministic  quantity  is  a  random  variable  which 
is  measurable  over  the  least  fine  b-field. 


*B  is  a  nonnull  atom  of  8  ,  if  PE  >  0 ,  and  B  contains  nc  other  sets  belong¬ 
ing  to  8  than  itself  and  the  empty  set. 


Tho  o- field  A  is  induced  by  the  atom  set  A  and  and  a  I iner  e-field 

6  containing  A  is  induced  by  the  atom  set  B^  ,  ,  B^  ,  B^  and  Bg 

in  Figure  28.  From  the  figure,  it  is  apparent  that  the  finer  the  o-field  8 
is,  the  closer  to  X(w)  the  conditional  expectation  E[X/B]  is.  If  a 
o-field  8  is  identical  to  the  o-field  induced  by  X(oj)  ,  then  i.tX/8]  is 
identical  to  X(w)  almost  surely.  The  variance  E[{X  -  E(X/8)}2]  is  pro¬ 
portional  to  the  area  between  two  random  variables  X(io)  and  £{X/8]  and 
the  conditional  variance  E [ { X  -  E[X/8]}2/8]  is  constant  on  every  nonnull 
atom  set  B  of  8  and  is  the  average  of  {X  -  E(X/6)}2  on  B  , 

Since  the  conditional  expectation  is  defined  in  an  equivalence 
sense,  the  area  under  the  conditional  expectation  E[X/8]  for  various  o~fi'ld 
8  must  be  identical  to  the  area  under  the  random  variable  X(w)  ,  which  i  ,  the 
unconditional  expectation  E[X]  .  The  conditional  expectation  E[X/8]  is 
the  closest  approximation  of  X(u)/  within  the  class  of  8-measurable  func¬ 
tions  in  the  sense  that  the  variance  is  minimized. 


APPENDIX  C 


RANDOM  WALK,  BROWNIAN  MOTION  AND  WHITE  NOISE 

Historically  random-walk  models  serve  as  a  first  approximation  to 
the  theory  of  diffusion  and  brownian  motion,  where  small  particles  are  ex¬ 
posed  to  a  tremendous  number  of  molecular  shocks.  Each  shock  has  a  negli¬ 
gible  effect,  but  the  superposition  of  many  small  actions  produces  an  obser¬ 
vable  motion.  Accordingly  we  want  to  present  a  random  walk  where  the  indi¬ 
vidual  steps  are  extremely  small  and  occur  in  very  rapid  succession.  In  the 
limit,  the  process  will  appear  as  a  continuous  motion,  i.e.,  the  so-called 
Brownian  motion.  Once  we  have  grasped  the  concept  of  Brownian  motion,  the 
white  noise,  which  is  fictitious  and  nonexsiting  but  enables  human  beings  to 
han  Le  many  mathematical  problems,  can  be  formally  defined  as  a  time  deriva¬ 
tive  of  Brownian  motion.  Here  a  brief  summary  of  Papoulis*  (36)  discussion 
on  the  subject  is  presented. 

The  underlying  experiment  is  the  tossing  of  a  fair  coin  an  infinite 
number  of  times,  and  each  tossing  occurs  every  T  seconds.  At  each  tossing 
we  take  a  step,  to  the  right  if  he  Is  show,  to  the  left  if  tails  show.  Our 
position  at  t  will  be  denoted  by  X(t)  .  Clearly,  X(t)  depends  on  the 
experimental  outcome,  i.e.  on  the  particular  sequence  of  heads  and  tails. 

We  have  thus  created  a  stochastic  process  known  as  random  walk.  Each  sample 
function  of  this  process  is  of  stair  ase  form  as  in  Fig.  29  with  discontin¬ 
uities  at  the  points  t  =  nT  the  steps  occur  instantly  and  their  length 
equals  S  . 

We  denote  by  x.  a  random  variable  equal,  to  .tS,  if  our  i-—'-  step 
is  to  the  right  or  lefr,  i.e.  heads  or  tails.  Thus 
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P{X.  = 
1 


=  p  =  h 


Plx.  = 

i 


=  q  =  h 


E{xi)  = 


E{ x? }  = 


Note  that  the  random  variable  x^  is  independent  and  has  zero  mean.  The 
position  at  t  =  nT  is  clearly  a  random  variable  given  by 


X(nT)  = 


X1  +  X2  + 


(C-l) 


Suppose  that  after  the  first  n  tossings,  k  heads  show,  then  the  value  of 
X(nT)  would  be  given  by 


XCnT)  =  kS  -  (n  -  k)S  =  (2k  -  n)S  =  rS 


(C-2) 


where 


r  =  2k  -  n 


Since  (X(nT)  =  rS}  is  the  event  {k 


r  +  n 


(C-3) 


heads  in  n  tossing}  the 


probability  is  given  by 


P{X(nT)  =  rS}  =  P heads}  = 


n  +  r 


(C-4) 


If  n  is  large  and  npq  »  1  ,  Demoivre -Laplace  theorem  (38}  is  applicable 
for  values  of  k  in  the  /npq  neighborhood  of  its  most  likely  value  np  , 


/npq  <  k  <  np  +  /npq 


and  the  approximate  form  of  Eq.  (C-4 )  is  given  by 


(C~5) 
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P{k  heads}  = 


1 _ .  _  (k  -  nP ); 


'll  n  'p  q 


2npq 


(C-6 ) 


Substitution  of  P  =  q  =  %  and  k 

following  approximate  expression 


n  +  r 


into  Eq.  (C-6)  yields  the 


P{X(nt)  =  rS}  2 


1  r* 

e  “  2n 


/n/2 


(C-7 ) 


provided  that  r  is  of  the  order  of  /n  . 

Furthermore ,  it  can  be  shown  that* 


P{X(nT)  <  rS}  =  h  +  erf 


(C-8) 


»n 


Finally,  the  mean  and  variance  of  the  random  variable  X(nT)  are  easily  ob¬ 
tained  and  they  are  as  follows: 


E[x(nT)]  =  0 

E(x2(nT)]  =  nS2 


(C-9 ) 
(C-10) 


In  the  following  discussion.  Brownian  motion  is  developed  as  a  .limiting  form 
of  the  random  walk.  For  the  time 


t  =  nT 


(C-ll) 


the  mean  and  variance  of  X(t)  become 


E[X(t )] 
E(X2(t>] 


=  0 


tsf_ 

"t 


(C-12) 
(C-13 ) 


*Athanasias  Papoulis,  Probability,  Random  Variables,  and  Stochastic  Processes, 
McGraw-Hill,  p.  08. 
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Suppose  now  that  we  keep  t  constant  but  we  make  C  and  T  tend  to  zero. 
The  variance  of  X(t)  will  remain  finite  and  different  from  zero  only  if  S 
tends  to  zero  as  /?  .  Otherwise,  X(t)  would  be  meaningless.  Therefore, 
assuming 

S2  =  aT  (C-14) 

wo  define  the  process  W(t)  as  a  limit 

W(t)  =  X(t)  (C-15) 

A  family  of  continuous  functions  results  for  almost  all  outcomes,  which  is 
known  as  a  Brownian  motion  or  Wiener-Levy  process.  From  5qs.  (C-9)  and 
(C-10 ) ,  the  mean  and  variance  of  this  process  are  obtained  and  they  are 

E[W(t)J  =  0  (C-16) 

E[W2 (t) ]  =  at  (C-17) 

The  value  of  random  process  W(t)  can  be  determined  from  dq.  (C-2)  and  given 
by 

W  =  rS  (C-18) 


In  connection  with  Eqs.  (C-J.l),  (C-14)  and  (C-18),  we  have  the  following  ex¬ 
pression: 


r_  _  W/S  _ _ W _  _  _W_ 

/n  /E7t  /tPTf  /at 


(C-19) 


and  hence  the  probability  distribution  F(>r,t)  is  obtained  as  a  limit  of 
Etj..(C-8) 


F(W.t)  =  ?{'W(r)  <  W,>  »  Vj  *  erf 


(C-20) 


/at 


The  probability  density  f(W.t)  is  readily  determined  from  Cqt  (C-20 )  arid 


is  given  by 


f(«tt)  =  ~z  GXP  (-  frr) 
/2at 


(C-21 ) 


Thus,  the  random  process  W(t)  is  normal,  with  zero  mean  and  variance  at  . 

The  fact  to  be  pointed  out  here  is  that  in  passing  to  a  limit,  all 
formulas  for  the  process  VJ(t)  remain  meaningful  and  agree  with  physically 
significant  formulas  of  diffusion  theory  which  can  be  derived  under  much 
more  general  conditions  by  more  streamlined  methods  (Einstein-rfiener  theory 
and  Uhlenbeck-Orstein  theory).  For  example,  the  density  funccion  is  obtained 
as  a  solution  of  the  diffusion  equation  by  Einstein  (37).  The  same  thing  can 
be  done  by  using  the  autocorrelation  of  the  solution  of  the  Langevine  equa¬ 
te  (36). 

From  Eq.  (C-l),  it  is  seen  that  for  t1  >  t2  »  W(t^)  “  W(t^)  is 

independent  of  W(t2)  -  W(0)  =  W(t^)  .  Hence, 

E[{W(t1)  -  W(t2)}W(t2)]  =  E(W(t1)  -  W(t2)]E[W(t2)]  =  0  (C-22) 

Thus, 

E[W(t1)W(t2)3  -  E[W2(t2)]  =  0  (C-23) 

Since  the  left  hand  side  of  Eq.  (C-22)  is  an  autocorrelation  R(t. ,t0)  and 

d.  z 

the  right  hand  side  is  ctt0  from  Eq.  (C-17),  the  following  is  developed. 


R(W  = 


ctt 2  for  t1  >  t2 


atj  for  tj  J  t2 


(C-24) 


An  infinitesimal  increment  d6t  of  Brownian  motion  W(t)  .is  defined  as 


d0t  =  W(t  +  e)  -  W(t) 


The  covariance  of  dBt  is  determined  from  Eq.  (C-24)  as  follows: 


(C-25) 
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Eld6t,dB  ]  =  l'.[{W(t  +  e)  -  W(t)}{W(t  +  e)  -  W(t)}] 

=  R(t  +  e,t  +  o)  -  R(t  e,t)  =  R(t,t  +  e)  +  R(  t ,  t ) 


=  a(t  +  e)  -  at  -  at  +  at 


E[dg  dgj  =  ae 
X  U 


(C-26) 


If  we  formally  define  white  gauss ian  process  u(t)  as  a  time  derivative  of 
Brownian  motion  as  follows : 


dW(t)  _  lim  Pt 
dt  e-K)  e 


(C-27) 


Then  the  autocorrelation  of  white  gaussian  process  u(t)  is  of  the  form 


E[u(t)u(t)]  =  =  a6(t  -  t) 


(C-29 ) 


Therefore,  the  variance  of  white  Gaussian  process  u(t)  is  infinite,  which 
is  in  agreement  with  the  axiomatic  definition. 


APPENDIX  D 


MARTINGALES 


if 


A  sequence  of  random  variables  xi*x2’’’*  ca^ec^  a  mar'tingale 


E[|xn|]  <  »  ,  n  > 


and 


( D— 1 ) 


E[x  l4/x, ,••• ,x  ]  =  x 

1  n+1  1  n  n 


with  Probability  1. 

A  stochastic  process  {x^jteT}  is  called  a  martingale  if 

E[|xt|]  <  » 


for  all  t  and  if,  whenever  n  >  1  and  t^  <  •••  <  t 


n+1 


E[x.  /x  ]  =  x 

n+1  tl  n  n 


with  Probability  1. 


(D-2) 


<Theorem  1> 


If  y1»y2»y3»***  are  defined  as 


y,  =  X,  ,  y,  =  x2  -  x!  »  ya  =  xa  "  x 


then,  if  the  x^  process  is  a  martingale. 


3  rt2  5 


(D-3) 


E E | yn 1 3  <  00  >  Etyn+i/yi’’‘',ynJ  =  0  »  n  -  1  »  (D-4) 

with  Probability  1.  The  x^’s  are  thus  partial  sums  of  the  series  l  y^  , 
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where  the  y^'s  satisfy  the  condition  (D-4).  Conversely,  the  partial  sums 


of  any  such  series  constitute  a  martingale. 


<Proof> 


xn’s  and  y^'s  are  linearly  related  and,  hence,  the  mapping  ma¬ 


trix  is  of  full  rank.  Therefore,  the  inverse  of  the  mapping  matrix  exists, 
which  implies  that  the  conditioning  {y^,***,y  ]  (D-4)  can  be  replaced 


by  the  conditioning  {x. ,***,x  }  . 

i  n 

Thus 

E[Wyl’‘"’yn'  = 

E[(x  , .  -  x  )/x, ,*•• ,x  ] 
ntl  n  1  nJ 

(D-5) 

E[Wyl’"->yn'  = 

E(x  ^./x. • ,x  ]  -  x 
n+1  1  n  n 

(D-6) 

Since  the  xn  process  is  a  martingale,  the  right  hand  side  of  £a.  (D-6)  be¬ 
comes  identically  zero  and  the  condition  (D-4)  immediately  follows.  Let 


x  =  y  y  and  the  y  's  satisfy  the  condition  (D-4),  then 
n  . n  n 

i=l 


"•*„<  = 

E(y1 

+  ••• +  Wxi>,”*xi.1 

(D-7) 

^Wv 

E[y1 

+  •"  *  Wyl’”’’yn’ 

(D— 8  ) 

E1WV 

*i  + 

•••  +  y  +  E[y  ^/y.  ,***,y  1 
lJn+l  Jn 

( D— 9  ) 

E'WV 

",xnl  = 

+ 

•  •  •  +  y  =  x 

Jn  n 

(D-10 ) 

,  the  x  '  s 

constitute  a  martingale  and  the  inverse  of  the 

theorem 

Let  n  ,5  ,  •  ♦ »  be  any  random  variables  with 


E[|n|]  <  « 


Then,  if  x  is  defined  by 
’  n  J 


x  =  E[n/C.  *•••  »5_] 
n  In 


(D-ll) 


the  process  is  a  martingale, 


< Proof > 


By  definition 


xn+l  * 


(B-12) 


Taking  the  conditional  expectation  of  Eq.  (D-12)  given  the  conditioning 
{£^,••*,5^}  ,  we  have  in  fact 

E[xn+l/5l’“‘,?n1  =  E[E{n/5l’"‘’5n+l}/Cl’“‘,5n]  (D"13) 

Since  the  o-field  generated  by  (C^,*** >^n+1)  contains  the  o-field  generated 
by  *  •  ,4^)  ,  the  smoothing  property  4  of  Appendix  B  can  be  applied  to 

the  right  hand  side  of  Eq.  (D-13),  Therefore,  Eq.  (D-13)  becomes 

E^Xn+l/?l’*"  =  »f’n^  =  xn  (D-14) 

with  Probability  1.  Since  are  random  variables  on  the  sample 

space  of  51»,,*t5n  » 

=  xn  (D‘15) 

with  Probability  1.  Taking  the  conditional  expectation  of  both  sides  of  Eq. 
(D-15)  given  {x^,*»*,xn}  and  using  the  smoothing  Property  4  of  Appendix  B, 
the  martingale  property  is  obtained. 


Therefore,  Theorem  2  is  proved. 


APPENDIX  r 


STOCHASTIC  FUNDAMENTAL  LEMMAS  AND  OFT IMALITY  CONDITION 


< Lemma  1> 


If  X  is  a  random  variable,  and  if 

E[XY]  =  0 


(F-l) 


for  every  deterministic  Y  (or  every  least  fine  o-field  measurable  random 
variable)  then, 


E[X]  =  0  a.s. 


(F-2) 


<Proof> 


E[XY]  =  0 
E[X]Y  =  0  a.s. 


since  Y  i  0 


E(X]  =  0  a.s. 


<Lemma  2> 


If  X  is  a  random  variable,  and  if 


E[XY]  =  0 


(F-3) 


for  every  to- function  Y  measurable  with  resepct  to  the  o-field  A  of  the 
measurable  w  set  ,  then 


E[X/A)  =  0  a.s. 


(F-4) 


< Proof > 

For  the  convenience  of  notation,  X  denotes  the  conditional  expec¬ 
tation  E[X/A]  . 
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APPENDIX  E 


TRACE  OF  MATRIX 

Let  A  be  an  (n*n)  square  matrix  and  a..^  represent  the  element 

of  it  1  row  and  jt^1  column,  then  the  trace  of  A  ,  denoted  by  tr(A)  ,  is 

T 

the  sum  of  the  diagonal  elements  of  A  ,  ana  similarily  of  A  ,  i.e., 

tr(A)  =  tr(A  )  =  a. „+•••+  a  =  j  a..  (E-l) 

11  .  nn  u 


'-Theorem  3> 


Let  A  be  any  (n*m)  matrix  and  represent  the  element  of 

,t.h  .  .th  T  T 

-i  '  row  ana  j  column,  then  AA  and  A  A  are  (nxn)  and  sa  lare 


matrices  respectively,  and  their  trace  is  uniquely  determined  by  the  sum  of 


square  of  elements,  a?.  ,  i.e., 


T  T  U  ™ 

tr(AA*)  =  tr(A  A)  =  £  la?. 

i«l  j=i  13 


(E-2) 


If  A  is  defined  to  be 


V  -  V 


(E-3) 


where  V  is  an  (n*m)  matrix  with  elements  and  V  is  an  approxima¬ 


tion  of  V  ,  with  elements  <f..  ,  then  A  represents  the  approximation  error 

3.  J 


with  elements  (v...  -  \h. )  and,  lienee,  tr  (AA*)  is  the  sum  of  the  square 


a  o  T 

errors,  (v. .  -  v..)z  .  Therefore,  tr(AA  )  would  be  a  sensible  criteria  tc 
13  3-3 


be  minimized  and  the  solution  V  is  the  least  square  error  solution.  If  V 


.T, 


is  a  matrix  of  random  variables  v..  and  the  risk,  trE[AA  1  is  minimized, 

J 


then  the  solution  V  is  the  minimum  variance  estimate.  Note  that  V  is  not 
necessarily  a  vector. 


E[XY]  =  E[(X  -  X  +  X)Y]  (F-5) 

Applying  the  smoothing  property  4,  Eq.  (F-5)  becomes 

E[XY]  =  E[E{(X  -  X)Y/A}]  +  E[XY/A]  a.s.  (F-6? 

E[XY]  =  E[XY]  =  0  a.s.  (F-7) 

Since  both  X  and  Y  are  A  measurable  random  variables  and  Y  is  arbi¬ 
trary,  it  is  possible  to  choose  Y  =  X  ,  then  Eq.  (F-7)  becomes 

E[XY1  =  E[XY]  =  E[X2]  =  0  a.s.  (F-8) 

The  above  is  true  only  for  X2  =  0  .  Therefore , 

X  =  E[X/A3  =  0  a.s. 


* 

£ 

1 

5 

I 

-A 
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i 

! 

1 
-  ^ 

i 

M 

(F-5) 

i 

M 

(f-6: 

1 

(F-7  ) 

1 

<Lemma  3> 


If  X(to,t)  is  a  stochastic  process  defined  on  the  set 


te[tlSt2]  ,  and  if 


E[X(u>,t)Y(u>,t)]dt  =  0 


(F-9) 


for  every  random  process  Y(oj,t)  measurable  with  respect  to  a  o-field  A(t)  , 


te[t1,t2]  »  then 


E[X(w,t)/A(t)3  =  0  a.s. 


(F-10) 


for  every  t  ,  tsJt^^jtj]  . 


<Proof> 


(*2 

j  F.lXYldt  =  0 

'  t 

1 


n^iWfWBWWfc 
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f  2  E[E{XV/A(t)}]dt 
\ 


=  0 


(F-ll) 


fX7 


E[E{X/A(t)}  •  Yjck  =  0 


(F-12) 


Suppose  that  E[X(u),t)/A(t)]  i-  0  .  Since  Y  is  arbitrary  A(t)  measur¬ 
able  fu  .  ion,  we  can  choose  Y  such  that 


Y(ui,t)  =  E[X(w,t)/A(t)] 


(F-13) 


Therefore,  Eq.  (F-12)  becomes 


f  2 

1  E[{E(X(a!,t)/A(t  )}2)dt  =  0  a.s. 


(F-14) 


The  above  is  positive  unless  E[X(w,t)/A(t)]  =  0  for  any  t  ,  teft^t^  . 

Therefore , 


E[X(ui,t)/A(t)]  =  0  a.s.  tett^jtj] 


(F-15) 


Let’s  consider  the  following  risk  function 

R(g)  =  frE((X  -  g)(X  -  g)T]  =  E[(X  -  g)T(X  -  g)J  (F-16) 


where  X  is  an  n*l  vector..  If  we  want  to  minimize  the  risk  (F-16)  with 
any  5-measurable  function  g  ,  then  the  solution  g  is  the  minimum  variance 
estimate.  In  order  to  minimize  .R  ,  we  introduce  a  variation  <$g  on  g  ,  i.e. 


g  =  g  +  <5g 


(F-17) 


Then  the  risk  (F-16)  becomes 

K(g  +  «g>  =  E[(X  -  g  +  dg)T(X  -  g  +  <5g)]  (F-18) 

Expanding  Eq.  (F-18)  about  g  ,  the  following  first  variation  is  obtained. 


18? 


«R(g)  =  E[(X  -  g)T6g]  +  E[6gT(X  -  g)} 


(F-19 ) 


For  the  minimum  of  R  ,  the  first  variation  6R  must  be  zero  and,  hence,  the 
following  must  be  satisfied 


E[(X  -  g)T<5g]  =  0 


( F- v0  > 


Since  a  variation  Sg  is  an  arbitrary  6-measurable  function,  the  stochastic 
fundamental  lemma  2  can  be  applied  and  the  optimality  condition  therefore  is 
obtained  as  follows: 


E[(X  -  g)/8}  =  0  a.s. 


(F-21 ) 


g  =  E[X/8] 


(F-22) 


The  solution  g  is  the  minimum  variance  estimate  and  is  given  as  the  condi¬ 
tional  expectation  of  X  given  B  . 

If  we  want  to  minimize  the  risk  (F— 16 )  with  a  linear  function  6  r 
Y  which  is  an  (mxl)  random  variable  observed,  i.e. 


g  =  KY 


(F-23) 


instead  of  any  S-measurable  function,  then  the  solution  g  is  the  linear 
minimum  variance  estimate.  This  time  we  have  to  determine  an  (nxm)  matrix 
K  within  a  class  of  B-measurable  functions  such  that  the  risk  (F-16)  is 
minimized.  Introducing  a  variation  6K  on  K  ,  the  risk  (F-16)  becomes 


R(K  +  <5K)  =  trE({X  -  (K  +  SK)Y}{X  -  <K  t  6K)Yr] 


(F-24) 


Expanding  Eq.  (F-24)  about  K  ,  the  first  variation  follows,  i.e. 

6R(K)  =  trE[(X  -  RY)YT6KT]  +  trE[6KY(X  -  RY)T}  =  0 


(F-25) 


Since  6K  is  an  arbitrary  B-measurable  function,  the  stochastic  fundamental 
Lemma  2  is  applied  and  the  optimality  condition  is  obtained. 
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E[(X  -  KY)YT/8]  =  0  a.s. 

(F-26) 

C[XYT/85  =  KE{YYT/B]  a.s. 

(F-27) 

K  =  E[XYT/8]{E[YYT/S] a.s. 

(F-28) 

Therefore,  the  linear  minimum  variance  estimate  is  given  by 

g  =  ecxy^shecyy^pd^y 

( F-29 ) 

For  the 

scalar  random  variables  X  and  Y  ,  the  linear  minimum 

variance  es- 

timate 

X  is  obtained  from  Eq.  (F-29 )  as  follows: 

%  =  asm, 

S[YY/83 

(F-30) 

If  we  choose  K  with  a  deterministic  number  which  is  measurable  over  the 
least  fine  a-field,  the  conditioning  becomes  unconditional  and  the  linear 
minimum  variance  estimate  (F-30)  becomes 
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